R/setup_muscles.r

In ime-luebeck/semgsim: Simulating surface electromyographic (sEMG) and force signals

#' Initialize the set of muscles to be simulated
#'
#' Convert data given in the config file to a data.frame.
#' 
#' @param muscles A list of muscle objects.
#' @return A data.frame containing a column 'muscle' giving the indices of the
#'   different muscles and a column 'muscle.obj' containing the corresponding
#'   muscle objects.
convert_muscle_list_to_df <- function(muscles) {

   stopifnot(is.list(muscles))
    
    for (muscle in muscles)
        stopifnot(class(muscle) == "muscle")

   data.frame(muscle = seq_along(muscles), muscle.obj = I(muscles))
}


#' Generate motor units (MUs) and muscle fibers for all muscles
#'
#' Set up the muscle structure by generating motor units (MUs) and their
#' contained muscle fibers according to the given configuration.
#'
#' Generates a set of MUs parameters following an exponential distribution with
#' muscle-specific parameters and then gathers all the MUs generated for the
#' different muscles in a single data.frame.
#' 
#' @param muscles A data.frame containing information on the muscle properties.
#'   Will usually be generated by a call to \link{convert_muscle_list_to_df}.
#'
#' @return A data.frame containing categorical columns 'muscle' and 'MU'
#'   giving the corresponding indices and 'MU.obj' containing the actual MU
#'   objects. These in turn have a property $fibers that contains a list of the
#'   muscle fibers contained in that MU. No geometrical information is present
#'   at this stage.
#' 
generate_MUs <- function(muscles) {

    stopifnot(is.data.frame(muscles))
    
    MUs_list <- lapply(seq_len(nrow(muscles)), function(i) {

        ## Generate MUs contained in the muscle, i.e., assign parameter values
        MUs <- generate_MUs_for_muscle(muscles$muscle.obj[[i]])

        ## Distribute MUs according to muscle geometry
        MUs <- distribute_MUs_in_muscle(muscles$muscle.obj[[i]], MUs)

        ## Convert to data.frame
        MUs <- data.frame(MU = seq_along(MUs), MU.obj = I(MUs))

        ## Populate the MUs with fibers
        for (j in seq_along(MUs$MU.obj))
            MUs$MU.obj[[j]]$fibers <-
                generate_fibers(muscles$muscle.obj[[i]], MUs$MU.obj[[j]])
            
        cbind(muscle = i, MUs)
    })

    MUs_df <- plyr::rbind.fill(MUs_list)
    stopifnot(is.data.frame(MUs_df))
    MUs_df
}


#' Generate Motor Units for a Muscle
#'
#' Assign MU properties according to muscle parameters.
#'
#' MUs have a set of parameters (recruitment thresholds, peak twitch force,
#' number of innervated fibers) that are closely related, and which follow an
#' exponential distribution: There are many MUs that are generated early, have a
#' low peak twitch force and innervate few fibers. On the other side of the
#' spectrum, there are few MUs that are recruited only at high force demand
#' levels, attain a high peak twitch force and innervate many fibers.
#' The exact parameters of the distribution are a property of the muscle.
#'
#' 
#' @param muscle A muscle object.
#' @return A list of \code{muscle$num_MUs} MU objects. Contains MU objects of
#'   class either "MU_rect" or "MU_circ", as defined by \code{muscle$MU_class}.
#' 
generate_MUs_for_muscle <- function(muscle) {

    stopifnot(
        class(muscle) == "muscle",
        muscle$shape == "rect"
    )

    ## Calculate general MU pool parameters
    rec_start <- muscle$MU_size_params$rec_thresh$min
    rec_stop <- muscle$MU_size_params$rec_thresh$max
    b <- muscle$rate_coding_params$shape_b
    a1 <- log(rec_stop) / muscle$num_MUs
    a2 <- log(rec_stop/ b) / muscle$num_MUs

    muscle$MU_size_params$rec_thresh$min <- rec_start
    
    MUs <- list()
    
    for (i in seq(1, muscle$num_MUs)) {

        MU <- list()
        class(MU) <- muscle$MU_class

        if (muscle$rate_coding_params$rec_model == 1)
            MU$rec_thresh <- exp(a1 * i)
        else
            MU$rec_thresh <- b * i * exp(a2 * i) / muscle$num_MUs
        
        ## Size principle parameters are monotonously related to the recruitment threshold
    		## size_factor is between 0 and 1 and scales properties between the smallest and
    		## the largest MU in the muscle
        size_factor <- (MU$rec_thresh - rec_start) / (rec_stop - rec_start)
        MU <- update_MU_size_params(MU, muscle$MU_size_params, size_factor)
    		MU$num_fibers <- round(MU$peak_force / MU$peak_fiber_force)
    		## Nominal number of fibers to be expected judging from MU size, exluding
    		## stochastic variations (this will be required later on if the MU needs to be resized
    		## because its size violates the muscle boundaries).
    		MU$num_fibers_nom <- MU$peak_force / MU$peak_fiber_force_nom
		
        ## All other parameters are assumed to be randomly distributed following
        ## a Weibull distribution
        for (attribute_name in names(muscle$random_MU_params))
            MU[[attribute_name]] <-
                sample_bounded(
                    rweibull,
                    list(n=1,
                         shape=muscle$random_MU_params[[attribute_name]]$shape,
                         scale=muscle$random_MU_params[[attribute_name]]$scale),
                    list(min=muscle$random_MU_params[[attribute_name]]$min,
                         max=muscle$random_MU_params[[attribute_name]]$max))

        MU_twitch <- calc_MU_force_twitch(MU)
        MU$total_twitch_power <- MU_twitch$total_twitch_power
        MU$force_twitch <- MU_twitch$force_twitch_fun
        
        MU$gain_fac_fun <- function(ISI) {
          force(MU$twitch_t_rise)
          normalized_fr <- MU$twitch_t_rise / ISI
          if (normalized_fr <= 0.4)
            gainfac <- 1
          else {
            r <- muscle$twitch_r
            c <- muscle$twitch_c
            gainfac <- 0.4 / normalized_fr / (1 - r) *
              (1 - r * exp((0.4 - normalized_fr) / c))
          }
          gainfac
        }
        
        MUs[[i]] <- MU
    }

    stopifnot(is.list(MUs))
    MUs
}


calc_MU_force_twitch <- function(MU) {
  ## Setup force twitch model; see doc for details
  thr <- MU$twitch_t_half_rel
  tri <- MU$twitch_t_rise
  tl <- MU$twitch_t_lead
  
  ## These should all be far below 1s.
  stopifnot(thr < 1 && tri < 1 && tl < 1)
  
  k <- log(2) / (thr - tri * log(tri + thr/tri))
  m <- k * tri
  p <- MU$peak_force * exp(-k * tri * (log(tri) - 1))
  
  twitch <- list()
  
  twitch$total_twitch_power <- p * k^(-m-1) * gamma(m+1)
  stopifnot(twitch$total_twitch_power > 0)
  
  twitch$force_twitch_fun <- function(t) (t >= tl) * p * t^m * exp(-k * (t-tl))
  
  twitch
}


update_MU_size_params <- function(MU, MU_size_params, size_factor, skip_param="rec_thresh") {
  stopifnot(size_factor <= 1)
  ## Size principle parameters are monotonously related to the recruitment threshold
  for (attribute_name in names(MU_size_params)) {
    
    if (attribute_name == skip_param)
      next
    
    ## Parameters defining the relationships determining this MU attribute
    attribute_params <- MU_size_params[[attribute_name]]
    
    ## Shape of the relationship between this parameter and the recruitment threshold
    if (attribute_params$relation == 1)
      relation <- function(x) x
    else if (attribute_params$relation == 2 && attribute_name == 'peak_fiber_force') {
      ## automatically determine the relative importance of the squared and linear terms, 
      ## such that the number of fibers as a function of recruitment threshold has a 
      ## positive slope over the whole range. (See paper notes for calculations.)
      ## The maximum should be attained at size_factor=2
      Fmin <- MU_size_params$peak_force$min
      Fmax <- MU_size_params$peak_force$max
      fmin <- attribute_params$min
      fmax <- attribute_params$max
      dF <- Fmax - Fmin
      df <- fmax - fmin
      x <- 2
      alpha <- peak_fiber_force_alpha(Fmin, dF, fmin, df)
      ## Note that this shifts the minimum of the parabola to -alpha/(2-2alpha).
      ## Left of the minimum, we extrapolate the minimum value constantly
      ## (to avoid unphysiological increases of fiber force for very low fiber numbers)
      relation <- function(x) (x >= -alpha / (2-2*alpha)) * (alpha*x + (1-alpha) * x^2)
    } else if (attribute_params$relation == 2)
      relation <- function(x) (x >= 0) * x^2
    else
      stop('Unknown size principle relationship provided.')
    
    ## Calculate the nominal value of this parameter based on the size principle
    attribute_val_proposal <- 
      calc_affine_param(min_val = attribute_params$min,
                        max_val = attribute_params$max,
                        ratio = relation(size_factor))
    
    ## If the parameter is modelled as a random variable, sample from a Weibull
    ## distribution centered at the computed nominal value
    if ('shape' %in% names(attribute_params) && attribute_val_proposal > 0) {
      MU[[paste(attribute_name, 'nom', sep='_')]] <- attribute_val_proposal
      attribute_val_proposal <- rweibull(n=1, shape=attribute_params$shape,
                                         scale = attribute_val_proposal)
    }
    if (attribute_params$round == TRUE)
      attribute_val_proposal <- attribute_val_proposal %>% round
    
    MU[[attribute_name]] <- max(attribute_params$lower_bound, attribute_val_proposal)
  }
  
  MU
}


peak_fiber_force_alpha <- function(Fmin, dF, fmin, df) {
	## If the following condition is not met, we run into problems because deal with muscle
	## boundaries by adjusting the corresponding MU's fiber count and then appropriately scaling
	## all other properties to such smaller MU. If all MUs have roughly the same number of fibers 
	## anyway, this does not really make any sense.
	nfibmin_nom <- Fmin / fmin
	nfibmax_nom <- (Fmin+dF) / (fmin+df)
	if (nfibmax_nom / nfibmin_nom < 10)
		stop("Expecting at least ten-fold fiber number increase from the smallest to the largest MU. Adjust ranges for peak_force and peak_fiber_force params in muscle config.")
	x <- 2
	alpha <- (Fmin / dF - fmin / df) / ((2*x-1)*Fmin/dF + x^2) + 1
	alpha
}
	
#' Calculate a variable that is affine in another variable, given that
#' variable's value.
#'
#' Calculates
#' `var2 = var2min + (var1i - var1min) / (var1max - var1min) * (var2max - var2min)`
#' where the quotient `(var1i - var1min) / (var1max - var1min)` corresponds to
#' function parameter `ratio`.
#' 
#' @param min_val Minimum value of the variable to calculate
#' @param max_val Maximum value of the variable to calculate
#' @param ratio Value between 0 and 1, indicating position between min and max
#' @return (Scalar) value of the affine dependent parameter
#'
calc_affine_param <- function(min_val, max_val, ratio) {
    stopifnot(list(min_val, max_val, ratio) %>% plyr::laply(is.nice.scalar) %>%
              all)
    min_val + ratio * (max_val - min_val)
}


#' Place MUs somewhere in the muscle area.
#'
#' Distribute MUs throughout the base muscle's assigned area, taking care of
#' assumed fiber densities, muscle boundaries, etc. This is a purely technical
#' function; the algorithm used for MU placement must be provided externally via
#' `muscle$MU_distributor`.
#' 
#' Returned objects have one very important property: $transform, which
#' transforms coordinates relative to the MU base region (a (0,1)x(0,1) square)
#' into coordinates relative to the muscle's base region (also a [0,1]x[0,1]
#' square). It is a function that takes 2D vectors with values between 0 and 1
#' in both dimensions. The same is true for its return values.
#'
#' @param muscle A single muscle object, must have element `MU_distributor`.
#' @param MUs A list of MU objects.
#' @return A list of modified MU objects with assigned geometrical properties.
#' 
distribute_MUs_in_muscle <- function(muscle, MUs) {
    stopifnot(is.list(MUs))
    stopifnot(is.list(muscle))

    ## Assign a geometrical location to the MUs
    MUs %<>% (muscle$MU_distributor)

    ## Assign transform function to each MU
    for (i in seq_along(MUs)) {
        MUs[[i]]$transform <- local({
            i <- i
            function(pos) {
                coord_x_rel <- pos[1] * MUs[[i]]$length_x
                coord_y_rel <- pos[2] * MUs[[i]]$length_y
                coord_x_abs <-
                    coord_x_rel + MUs[[i]]$pos[1] - MUs[[i]]$length_x / 2
                coord_y_abs <-
                    coord_y_rel + MUs[[i]]$pos[2] - MUs[[i]]$length_y / 2
                c(coord_x_abs, coord_y_abs)
            }
        })

        if (MUs[[i]]$pos[1] - MUs[[i]]$length_x / 2 < 0 ||
            MUs[[i]]$pos[1] + MUs[[i]]$length_x / 2 > 1 ||
            MUs[[i]]$pos[2] - MUs[[i]]$length_y / 2 < 0 ||
            MUs[[i]]$pos[2] + MUs[[i]]$length_y / 2 > 1) {

            ## Some part of the MU region is cutoff because it lies in the
            ## vicinity of the muscle boundary. Accordingly, we reduce the 
            ## number of fibers in order to not increase fiber density near the
            ## muscle boundary.
            MUs[[i]]$cutoff_ratio <- MUs[[i]] %>% calc_MU_cutoff_ratio
			## This is the _nominal_ number of fibers that would be expected irrespective
			## of stochastic variability. The exact number of fibers is calculated further below.
            MUs[[i]]$num_fibers_nom <-
                ((1 - MUs[[i]]$cutoff_ratio) * MUs[[i]]$num_fibers_nom)
            
			## This is now a smaller MU, so we want to adjust all other MU properties such 
			## that they fit such a smaller MU! To do this, we recalculate the size_factor
			## such that it fits with the nominal number of fibers, cf. above.
            ## size_factor can be negative if the adjusted number of fibers is
            ## less than muscle$num_fibers_min. This is a feature, not a bug!
			## To calculate the size_factor, we solve 
			##    num_fibers = peak_force(size_factor) / peak_fiber_force(size_factor)
			## for size_factor.
			## peak_force(..) is a linear function, but peak_fiber_force(..) is usually
			## quadratic, hence we need to solve a quadrativ equation here.
			## See paper notes for the details of this derivation.
			peak_force_min <- muscle$MU_size_params$peak_force$min
			peak_force_max <- muscle$MU_size_params$peak_force$max
			delta_force <- peak_force_max - peak_force_min
			peak_fiber_force_min <- muscle$MU_size_params$peak_fiber_force$min
			peak_fiber_force_max <- muscle$MU_size_params$peak_fiber_force$max
			delta_fiber_force <- peak_fiber_force_max - peak_fiber_force_min
			num_fibers_nom <- MUs[[i]]$num_fibers_nom
			
			alpha <- peak_fiber_force_alpha(peak_force_min, delta_force, peak_fiber_force_min, delta_fiber_force)
			
			if (muscle$MU_size_params$peak_fiber_force$relation == 1 || abs(alpha - 1) < 1e-2)
			  ## Linear model (old, obsolete?)
			  ## If the second condition is true, we're *almost* linear; the quadratic inversion becomes ill-conditioned.
			  ## Thus switch to the linear derivation in this case.
			  size_factor <- (peak_force_min - num_fibers_nom * peak_fiber_force_min) / 
				(num_fibers_nom * delta_fiber_force - delta_force)
			
			else if (muscle$MU_size_params$peak_fiber_force$relation == 2) {
			  ## Quadratic model	
			  radicand <- (alpha * num_fibers_nom * delta_fiber_force - delta_force)^2 -
				4 * (1-alpha) * num_fibers_nom * delta_fiber_force * (num_fibers_nom * peak_fiber_force_min - peak_force_min)
			  if (radicand < 0)
				sqrt_term <- 0
			  else
				sqrt_term <- sqrt(radicand)
			  size_factor <- (delta_force - alpha * num_fibers_nom * delta_fiber_force - 
								sqrt_term) / (2 * (1-alpha) * num_fibers_nom * delta_fiber_force)
			} else
			  stop('Invalid size principle relation specified')
			  
			## Prop_factors > 1 yield very unpredictable and often unphysiological results.
			size_factor <- min(1, size_factor)
            MUs[[i]] <- update_MU_size_params(MUs[[i]], muscle$MU_size_params, size_factor, 
                                              skip_param = "")
								
			## If things went well, we should now have gotten nominal (average) values for
			## peak_force and peak_fiber_force that yield exactly num_fibers_nom.
			## (That was the whole reason for all the calcuations above.)
			stopifnot(abs(MUs[[i]]$peak_force / MUs[[i]]$peak_fiber_force_nom - MUs[[i]]$num_fibers_nom) < 1 || size_factor < 0)
								
			## The above call to update_MU_size_params has now re-sampled peak_force and
			## peak_fiber_force using the new size_factor. Using these two, we can now calculate
			## the actual fiber count, which may differ from the nominal fiber number due to random
			## variations in the two peak force quantities.
			## This is the actual number of fibers that we will generated later on.
			## (See distribute_fibers_unif.MU_circ.)
            MUs[[i]]$num_fibers <- round(MUs[[i]]$peak_force / MUs[[i]]$peak_fiber_force)
            
            ## Update force twitch properties accordingly
            MU_twitch <- calc_MU_force_twitch(MUs[[i]])
            MUs[[i]]$total_twitch_power <- MU_twitch$total_twitch_power
            MUs[[i]]$force_twitch <- MU_twitch$force_twitch_fun            
            
            stopifnot(list(MUs[[i]]$cutoff_ratio,
                           MUs[[i]]$num_fibers,
                           MUs[[i]]$peak_force,
                           MUs[[i]]$rec_thresh) %>%
                      lapply(is.nice.scalar) %>%
                      simplify2array %>%
                      all)
            stopifnot(MUs[[i]]$cutoff_ratio >= 0)
            stopifnot(MUs[[i]]$cutoff_ratio < 1)
            stopifnot(MUs[[i]]$num_fibers > 0)
            stopifnot(MUs[[i]]$peak_force > 0)
            stopifnot(MUs[[i]]$rec_thresh > 0)
        } else
            MUs[[i]]$cutoff_ratio <- 0
    }

    stopifnot(is.list(MUs))
    ## Reorder the list, since MUs may no longer be in ascending order of power due to
    ## adjustments during geometrical arrangements (see above).
    MUs[order(plyr::laply(MUs, function(MU) MU$rec_thresh))]
}


#' Distribute motor units on a regular grid across the muscle cross-section
#'
#' Deterministically create a regular grid of MU positions across the muscle
#' cross section for a muscle with a rectangular base shape. Assign the same
#' area
#' 
#' The MU areas are obtained by simple division of the unit square into a
#' (num_MUs_x x num_MUs_y) grid. No randomness is involved in the MU
#' distribution. All MUs are of the same size. Currently, only exact division of
#' the muscle area into MU areas is implemeted, i.e., no overlapping MUs and no
#' muscle sections without MUs.
#'
#' @param MUs A list of MU objects. Must contain \code{num_MUs_x * num_MUs_y}
#'   MUs.
#' @param num_MUs_x The number of MUs to be placed in one row in x direction.
#' @param num_MUs_y The number of MUs to be placed in one row in y direction.
#' @param overlap The amount of overlap to occur between neighboring MU regions.
#'   Currently, only \code{overlap == 0} is implemented.
#' @return A list of MU objects, containing the elements of \code{MUs} but with
#'   additional properties:
#' \itemize{
#'  \item{"$pos"}{The relative 2-D position (i.e. from [0,1] x [0,1]) of the MUs
#'                center on the muscle cross section.}
#'  \item{"$length_x"}{The relative size of the MU's surrounding rectangle box in
#'                   x direction on the muscle cross section. Value between 0
#'                   and 1. Used to transform coordinates relative to the MU
#'                   area to coordinates relative to the muscle area.}
#'  \item{"$length_y"}{The relative size of the MU's surrounding rectangle box in
#'                   y direction on the muscle cross section. Value between 0
#'                   and 1. Used to transform coordinates relative to the MU
#'                   area to coordinates relative to the muscle area.}
#' }
#' 
distribute_MUs_reg <- function(MUs, num_MUs_x, num_MUs_y, overlap = 0) {

    stopifnot(
        is.list(MUs),
        length(MUs) == num_MUs_x * num_MUs_y,
        class(MUs[[1]]) == "MU_rect" || class(MUs[[1]]) == "MU_circ",
        is.nice.scalar(num_MUs_x),
        is.nice.scalar(num_MUs_y),
        is.nice.scalar(overlap)
    )
    
    if (overlap != 0)
        fail("Overlapping of regular grids of MUs is not (yet) implemented.")
    
    ## Assign center positions on the muscle to all MUs.
    ## Positions are taken on a regular grid.

    step_x <- 1 / num_MUs_x
    step_y <- 1 / num_MUs_y

    ## This yields a (num_MUs_x * num_MUs_y) x 2 numeric containing 2-D tile
    ## center coordinates. 
    tile_center_coords <-
        expand_grid_vec(seq(from = step_x / 2,
                            to = 1 - step_x / 2,
                            by = step_x),
                        seq(from = step_y / 2,
                            to = 1 - step_y / 2,
                            by = step_y))

    for (i in seq_along(MUs)) {
        MUs[[i]]$pos <- tile_center_coords[i, ]
        MUs[[i]]$length_x <- step_x
        MUs[[i]]$length_y <- step_y
        MUs[[i]]$cutoff_ratio <- 0
    }
    
    MUs
}


distribute_MUs_unif <- function(MUs, muscle, tiling = c(1,1)) {

    stopifnot(
        is.list(MUs),
        class(MUs[[1]]) == "MU_rect" || class(MUs[[1]]) == "MU_circ",
        is.vector(tiling),
        length(tiling) == 2,
        is.nice.vector(tiling),
        all(floor(tiling) == tiling)
    )

    step_x <- 1 / tiling[1]
    step_y <- 1 / tiling[2]

    tile_corner_coords <- expand_grid_vec(seq(from = 0,
                                              by = step_x,
                                              to = 1 - step_x),
                                          seq(from = 0,
                                              by = step_y,
                                              to = 1 - step_y))
    
    for (i in seq_along(MUs)) {

        tile <- i %% (tiling[1] * tiling[2]) + 1

        MUs[[i]]$pos <- c(runif(n = 1,
                                min = tile_corner_coords[tile, 1],
                                max = tile_corner_coords[tile, 1] + step_x),
                          runif(n = 1,
                                min = tile_corner_coords[tile, 2],
                                max = tile_corner_coords[tile, 2] + step_y))
        
        MU_area <- MUs[[i]]$num_fibers / MUs[[i]]$fiber_density
        lengths <- calc_MU_lengths(MUs[[i]], MU_area, muscle)
        MUs[[i]]$length_x <- lengths$length_x
        MUs[[i]]$length_y <- lengths$length_y
    }

    MUs
}


#' Generate muscle fibers for a motor unit
#'
#' Take a muscle and a MU configuration and return a list of associated muscle
#' fibers for that MU.
#'
#' Initializes the geometry by calling \code{muscle$rel_fiber_distributor(MU)}
#' (which should be defined in the config file) and then calculates all other
#' fiber properties from the MU's properties and the geometrical position of the
#' fiber.
#'
#' @param muscle A muscle object.
#' @param MU A MU object.
#' 
#' @return A list of muscle fiber objects generated by applying the definitions
#'   given in the muscle and MU arguments.
#'
generate_fibers <- function(muscle, MU) {

    ## Place the fibers in the MU base region
    fiber_positions <- MU %>% (muscle$rel_fiber_distributor)

    stopifnot(is.list(fiber_positions),
              is.list(fiber_positions[[1]]),
              is.vector(fiber_positions[[1]]$pos_rel_to_MU),
              is.vector(fiber_positions[[1]]$pos_rel_to_muscle))

    ## Create fiber objects containing these positions.
    fibers <- fiber_positions %>% seq_along %>%
        lapply(function(id) {
            fiber <- list()
            class(fiber) <- 'fiber'
            fiber$base_pos_rel_to_muscle <-
                fiber_positions[[id]]$pos_rel_to_muscle

            ## Calculate the absolute positions of the fiber ends. These are
            ## (still) not final though, as they will be subject to some random
            ## variations.
            fiber$base_pos_abs <-
                fiber$base_pos_rel_to_muscle %>% (muscle$shape_transformation)
            fiber
        })
    
    ## Determine the remaining fiber properties.
    fibers <- lapply(fibers, function(fiber) {

        ## Determine the rotation applied to this fiber.
        fiber$rotation <- muscle$fiber_rotations(fiber$base_pos_rel_to_muscle)

        fiber$shape <- muscle$fiber_shape
        
        fiber$length_base <-
            (fiber$base_pos_rel_to_muscle) %>%
                (muscle$fiber_base_length_calculator)
        fiber$length <- fiber$length_base * muscle$fiber_length_randomizer()

		## Correct base position for randomness in length
		fiber$base_pos_abs <- fiber$base_pos_abs + (fiber$length_base - fiber$length) / 2 *
            (fiber$shape$pos(1) %>% rotate3D_lh(fiber$rotation))

        ## Determine the relative position of the fiber innervation zone.
        fiber$innervation_zone_rel <- muscle$innervation_zone_placer(fiber)

        ## Calculate the absolute position of the fiber innervation zone.
        fiber$innervation_zone <- fiber$base_pos_abs + fiber$length *
            (fiber$shape$pos(fiber$innervation_zone_rel) %>%
                 rotate3D_lh(fiber$rotation))

        fiber$end_point <- fiber$base_pos_abs + fiber$length *
            (fiber$shape$pos(1) %>% rotate3D_lh(fiber$rotation))
        
        ## Check that guarantees given in the config file hold
        stopifnot(fiber$length <= muscle$fiber_length_max,
                  max(fiber$innervation_zone_rel,
                      1 - fiber$innervation_zone_rel) <=
                          muscle$fiber_rel_half_length_max)
        
        fiber
    })

    fibers
}



## Take a MU object and return a list of 2-D positions relative (in the range
## (0,1) in both dimensions) to the MU base region.

distribute_fibers_reg <- function(MU)
    UseMethod("distribute_fibers_reg", MU)

distribute_fibers_reg.MU_rect <- function(MU) {

    ## Create a list of fiber positions relative to the (0,1)x(0,1) MU base
    ## square. Positions are taken from a regular grid.
    n <- ceiling(sqrt(MU$num_fibers))
    dt <- 1/(n+1)

    ## Here we create the actual coordinate grid
    coords <- expand_grid_vec(seq(dt, 1 - dt, dt), seq(dt, 1 - dt, dt))

    ## Split this numeric array into a list of its rows
    fiber_positions_rel_to_MU <- split(coords, row(coords))[1:MU$num_fibers]

    fiber_positions <- fiber_positions_rel_to_MU %>%
        lapply(function(pos_rel_to_MU) {
            pos_rel_to_muscle <- pos_rel_to_MU %>% (MU$transform)
            list(pos_rel_to_MU = pos_rel_to_MU,
                 pos_rel_to_muscle = pos_rel_to_muscle)
        })

    ## Return these relative positions.
    fiber_positions
}


distribute_fibers_reg.MU_circ <- function(MU)
    fail("distribute_fibers_reg.MU_circ not implemented.")

distribute_fibers_unif <- function(MU)
    UseMethod("distribute_fibers_unif", MU)

distribute_fibers_unif.MU_rect <- function(MU)
    fail("distribute_fibers_unif.MU_rect not implemented")

distribute_fibers_unif.MU_circ <- function(MU) {

    circ_midpoint <- c(0.5,0.5)

    fiber_count <- 0
    fiber_positions <- list()

    while(fiber_count < MU$num_fibers) {

        pos_rel_to_MU <- sample_points_in_circ_unif(1)

        pos_rel_to_muscle <- pos_rel_to_MU %>% (MU$transform)

        if (!(pos_rel_to_muscle %>% is_outside_unit_box)) {

            fiber_count <- fiber_count + 1
            fiber_positions[[fiber_count]] <-
                list(pos_rel_to_MU = pos_rel_to_MU,
                     pos_rel_to_muscle = pos_rel_to_muscle)
        }
    }

    fiber_positions
}

calc_MU_lengths <- function(MU, MU_area, muscle)
    UseMethod("calc_MU_lengths", MU)    

calc_MU_lengths.MU_circ <- function(MU, MU_area, muscle) {

    lengths <- list()
	
	# We want the ellipsis to have the same axis ratio as the whole muscle.
	# (E.g., if the muscle is flat and wide, we also want the MU area to be flat and wide.)
	# So we use a circle in normalized ([0,1]x[0,1]) muscle coordinates, which then
	# automatically becomes an ellipsis with the desired axes ratio in real/cartesian coordinates.
	# Radius r in normalized coordinates becomes
	# x_radius=r*muscle_length_x and y_radius=r*muscle_length_y
	# in real coordinates.
	# The area of an ellipsis is A=pi*x_radius*y_radius, hence we have
	# A=pi*x_radius*y_radius=pi*r^2*muscle_length_x*muscle_length_y
	# and thus we need the following radius in normalized coordinates:
	r = sqrt(MU_area/pi/muscle$width/muscle$thickness)
	lengths$length_x <- r
	lengths$length_y <- r
	
	# OLD, OBSOLETE (and possibly also wrong?)
    #radius <- sqrt(MU_area / pi)
    # 
    ### The rescaling by the muscle dimensions ensures a circular MU in the
    ### 'real' coordinate system.
    #lengths$length_x <- 2 * radius / muscle$width
    #lengths$length_y <- 2 * radius / muscle$thickness
    lengths
}

calc_MU_lengths.MU_rect <- function(MU, MU_area, muscle)
    fail("Not implemented")


calc_MU_cutoff_ratio <- function(MU)
    UseMethod("calc_MU_cutoff_ratio", MU)

calc_MU_cutoff_ratio.MU_circ <- function(MU, num_samples = 10000) {
    coords <- sample_points_in_circ_unif(num_samples)
    coords_transformed <- coords %>% plyr::aaply(.margins = 1,
                                                .fun = MU$transform)
    cutoff_mask <- is_outside_unit_box(coords_transformed)
    cutoff_ratio <- sum(cutoff_mask) / num_samples
    cutoff_ratio
}

calc_MU_cutoff_ratio.MU_rect <- function(MU)
    fail("Not implemented.")