# Force Feedback Control For Dexterous Robotic Hands Using Conditional
Postural Synergies

Dimitrios Dimou1 José Santos-Victor1 Plinio Moreno1 1Institute for Systems and
Robotics, Instituto Superior Tecnico, Universidade de Lisboa, Portugal.
Emails<EMAIL_ADDRESS>{jasv<EMAIL_ADDRESS>
###### Abstract

We present a force feedback controller for a dexterous robotic hand equipped
with force sensors on its fingertips. Our controller uses the conditional
postural synergies framework to generate the grasp postures, i.e. the finger
configuration of the robot, at each time step based on forces measured on the
robot’s fingertips. Using this framework we are able to control the hand
during different grasp types using only one variable, the grasp size, which we
define as the distance between the tip of the thumb and the index finger.
Instead of controlling the finger limbs independently, our controller
generates control signals for all the hand joints in a (low-dimensional)
shared space (i.e. synergy space). In addition, our approach is modular, which
allows to execute various types of precision grips, by changing the synergy
space according to the type of grasp. We show that our controller is able to
lift objects of various weights and materials, adjust the grasp configuration
during changes in the object’s weight, and perform object placements and
object handovers.

##  I Introduction

To perform complex manipulation tasks in unstructured environments, humans use
tactile feedback from their fingers. This feedback is provided by tactile
afferents located in the skin of the hand. Particularly, for handling small
objects with precise movements, the afferents located in the fingertips are
used, which have high density and adapt fast to pressure changes [1]. These
afferents provide information about the characteristics of the exerted contact
forces, such as the magnitude and the direction. For anthropomorphic robots to
be able to perform dexterous tasks similar force feedback signals must be used
to alleviate problems arising from uncertainty in measurements, and handle
external perturbations. For example, using open-loop position control to lift
a heavy object may fail due to slip without any feedback mechanism to provide
tactile information.

Previous works have used tactile sensors to design force controllers that use
slip prediction to update the desired normal forces applied by the fingertips.
The slip predictors are based on machine learning models such as neural
networks and random forests to classify multi-modal signals from a tactile
sensor. In all previous works, each finger was separately controlled by an
independent force controller. In addition, they required labeled data to train
the slip predictors and because each finger is controlled independently is not
obvious how to implement different anthropomorphic grasp types.

In this work we develop a force controller that takes as input the force
readings of the fingertips and computes the grasp size which is then used
along with a grasp type label to generate a grasp posture with the desired
characteristics. To avoid slippage the desired normal contact force is
calculated to be proportional to the tangential contact forces. The applied
normal force is then controlled using the size of the grasp as a control
variable. Larger grasp sizes mean less force is applied to the object. So the
grasp size is calculated from the error between the desired normal force and
the actual measured normal force. The grasp size is then given to the posture
sampler that generates a grasp posture, i.e. the finger joint angles. The
posture sampler is modeled with a conditional Variational Auto-Encoder (cVAE)
based on the framework proposed in [2]. With this framework we abstract away
the low-level control of the fingers and generate hand postures based on high-
level properties such as the type and the size of the grasp. So it works as a
mapping function that takes as input a low-dimensional vector and the grasp
type and size as conditional variables and maps them to a set of joint angles.

We show that with our controller we can control a dexterous robotic hand to
lift objects of different weights using three precision grasps. Our controller
is also able to compensate and retain a stable grasp during changes in the
objects’ weight, for example when filling up a cup or emptying it. In addition
we show how with the addition of the hand pose information we can use the
controller to calculate if the tangential force is due to gravity or due to a
support surface and use this information to perform handovers and place down
objects on surfaces. We perform several real-world experiments with a
dexterous robotic hand to showcase the capabilities of our controller and
support our design choices.

To sum up our main contributions are

  * •

We develop a controller for a dexterous robotic hand that uses force feedback
and the conditional synergies framework to perform dexterous manipulation
tasks.

  * •

We show that with our controller we can easily use different precision grasp
types, by changing only the grasp type variable which is given to the grasp
posture mapping function.

  * •

We demonstrate by incorporating information about the world pose of the hand
we can use our controller to perform additional tasks such as placing down and
handing over objects.

##  II Related Work

Roboticists have looked for inspiration in humans for developing methods for
complex object manipulation [3]. Neuroscientists have studied for a long time
the processes that allow humans to use tactile feedback to perform complex
manipulation tasks. Humans tend to adjust the grip force according to the
object’s weight, its friction and they use a safety margin to account for
uncertainties [4]. To gather information about the tactile states they use
multiple afferents that are located in the skin of the fingers [1]. There are
different afferents in different parts of the hand depending on their usage,
e.g. fast adapting afferents in the fingertips for precise manipulation. Based
on signals from these afferents, humans encode simple contact events into
action phases, such as grasping, lifting or releasing, which they combine in
order to perform more complex and long-horizon manipulation tasks [5].

In robotics tactile sensors have been used for object stabilization and slip
prediction in a variety of settings. For example, in [6], a compliant
anthropomorphic prosthetic hand was controlled using force sensing to maintain
object stability and avoid slip. In [7], they develop a control approach that
uses integrated force and spatial tactile signals to avoid slip with unknown
objects in real world settings. In [8, 9], grasp quality metrics are computed
based on the tactile feedback from the robots fingertips. In these works,
simple two or three fingered grippers were considered for simple grasping
tasks.

Force control with anthropomorphic robotic hands has also been explored in
more recent works. In [10], they employ three slip prediction methods to
estimate when slip starts and based on the force signals at that moment they
calculate the friction coefficient value. Based on the calculated friction
coefficient, they design a force controller that independently controls each
finger to achieve a desired normal force. The desired normal contact force is
set to be proportional to the tangential contact force and a safety margin
based on the evidence found in [4]. In [11], they train a random forest to
classify the contact states into the classes: no contact, contact, slip. Based
on this classification signal, when slip is detected they increase the desired
normal contact force to avoid it. In [12] they train a recurrent neural
network to estimate slip and the object material from the readings of a Biotac
sensor. The force controller is increasing the desired normal contact force
when slip is detected. All these works [12, 11, 10] use tactile feedback
sensors to predict slip. They collect labeled data, on which they train their
models. This approach is based on complex and expensive tactile sensors, and
the process of collecting data is cumbersome. In addition, the data do not
cover all possible hand poses, which would be impractical.

In contrast, in our work we do not rely on slip prediction, we avoid slip by
defining a tangential force gain and a safety margin that work for a large
number of objects. Furthermore, instead of independently controlling each
finger we use a synergistic framework to generate grasp postures, that is
conditioned on two variables: the grasp type and the grasp size. This way,
instead of controlling the values of each joint of each finger, we control
only the two conditional variables greatly simplifying the control pipeline.
This also, gives us the ability to use different grasp types in our
manipulation tasks by changing only the grasp type variable. In [13] also a
synergistic framework was used to prevent an object from slipping from a
humanoid hand, but they modeled only one synergy for a tripod grasp and they
used the forces on the robotic arm as feedback, while we use force feedback
from the fingertips. Our control algorithm could also be applied to different
hands as it does not depend on the hands configuration. Finally, in previous
approaches only lifting tasks had been considered. In our work we demonstrate
that our approach can be used to perform more complex tasks, such as placing
objects on surfaces and performing handovers, which was not done in previous
works.

Figure 1:  Example of modeling the contacts and friction during manipulation.

##  III Methods

Our goal in this work is to design a control algorithm for an anthropomorphic
robotic hand to perform dexterous manipulation skills such as lifting and
placing down objects. Our control algorithm will use tactile feedback from the
force sensors on the fingertips of the hand to decide the forces that need to
be applied to the object in each step of the task. Given the desired forces to
be applied, the size of the grasp will be computed. Given the grasp size and a
desired grasp type, the posture generator will generate a grasp posture, i.e.
the hand configuration, such that the force constraints are satisfied.

To model the contacts and friction we use Coulombs’ law, which states that in
order to avoid slip, the normal contact force $f_{n}$ to the contact surface
of an object, times the fiction coefficient $\mu$, has to be larger than the
tangential force $f_{t}$ [14]:

$\mu f_{n}\geq f_{t}$

You can see an example in Figure 1, where an object is pressed against a wall
by an applied normal force $f_{n}$, and we have the tangential force
$f_{t}=mg$ due to gravity. In order for the object to remain stable we need to
apply a normal force:

$f_{n}\geq\frac{f_{t}}{\mu}$

where $\mu$ is the friction coefficient between the object and the wall. In
the case of a dexterous hand manipulating an object, we want the normal forces
applied by all fingers to be greater than the tangential force divided by the
friction coefficient of the materials of the object and the fingertip.

Figure 2:  Schematic representation of the proposed force controller. The
input is the state (GRASP or RELEASE) and the force readings. Based on that
the grasp size is adjusted by a value $C$ and is given to the posture mapping
function along with the desired grasp type. A finger configuration is then
generated and commanded to the robot.

Since it is hard to accurately compute the friction coefficient between all
possible object materials previous works have used multi-modal tactile sensors
like the BioTac sensor, which provides information about the pressure, skin
deformation, and temperature, to predict slip and based on that signal to
increase the applied normal force. In our work we use the FTS3 sensors [15]
which is a low-cost sensor that measures the 3D force applied in each
fingertip. In addition, previous works gathered labeled datasets in order to
train their slip prediction models which is time-consuming and limits the
possible orientations of the hand, because gathering labeled data for all
possible orientations is impractical. To overcome this we experimentally
selected the parameters that determine the value of the applied normal force
such that we avoid slip for all objects in our dataset, from the lightest to
the heaviest.

In order to guarantee contact between the fingertip and the object, in the
beginning of the grasping phase, we use an offset $f_{n}^{offset}$ as the
minimum normal force applied by each finger. In [4] they also suggest that
humans use an additional safety margin which is proportional to the tangential
force, $f_{n}^{margin}\propto f_{t}$. So the final desired normal contact
force becomes:

$f_{n}^{des}=G\cdot f_{t}+f_{n}^{offset}$

where $G$ is the gain that includes the friction coefficient and the
additional safety margin.

To alleviate the effects of noise in the sensors, the running average of the
measured normal force $f_{n}$ and tangential force $f_{t}$ is used, as a low
pass filter. So for each force measurement we have the following relation:

$\hat{f}(n+1)=\alpha f(n+1)+(1-\alpha)f(n)$

where $\alpha\in(0,1)$ is a parameter that determines how much new
measurements affect the value, and is experimentally selected.

Figure 3:  Our control algorithm in Python-like pseudocode.

Given the measured normal force $f_{n}$ from the fingertip sensors we can
compute the error $f_{n}^{err}=f_{n}^{des}-f_{n}$. We use this error signal to
control the grasp size variable $g_{size}$, that we use as a conditional
variable in our posture mapping function. The grasp size represents the
distance between the thumb and the index finger in a grasp posture. So a
smaller grasp size will result in a tighter grasp and greater normal force
applied to the surface of the object. We use a linear controller for the grasp
size variable that is implemented as follows:

$g_{size}(n+1)=g_{size}(n)-K\cdot f_{n}^{err}$

where $K$ is a parameter that controls the rate of decrease of the grasp size,
and is experimentally selected. So when the error between the desired normal
force and the actual normal force is large the grasp size decreases so tighter
grasp postures are generated in order to apply more normal force. In practice,
in order to avoid oscillations in the grasp size we use the desired normal
force as a high threshold that we want the measured normal force to be below:

$f_{n}<f_{n}^{des}=f_{n}^{threshold\\_high}.$

If the normal force is below that threshold the grasp size does not change
even if there are small oscillations in the measured tangential and normal
forces. Also, in order to avoid the hand applying too much force that damages
the hardware or the object we use a low threshold, that is:

$f_{n}>f_{n}^{threshold\\_low}=f_{n}^{threshold\\_high}-w_{threshold},$

where $w_{threshold}$ is the width of the threshold in $mN$. If the measured
normal force is below the grasp size increases in order to apply less force.
So the final grasp size variable for grasping is calculated as follows:

$g_{size}(n+1)=\begin{cases}g_{size}(n)-K\cdot f_{n}^{err\\_high}&\text{if
$f_{n}^{err\\_high}>0$}\\\ g_{size}(n)+K\cdot f_{n}^{err\\_low}&\text{if
$f_{n}^{err\\_low}<0$}\end{cases}$ (1)

$\text{where}\begin{cases}f_{n}^{err\\_high}=f_{n}^{threshold\\_high}-f_{n}\\\
f_{n}^{err\\_low}=f_{n}^{threshold\\_low}-f_{n}\end{cases}$

This is similar to the deadband control method [16], where instead of having a
fixed reference point, an operating range is set. If the response is in this
range, the controller does not exert any correction. In our case, the
operating range changes according to the force signals from the robot’s
fingertips.

The grasp posture mapping function is based on the conditional postural
synergies model presented in [2]. It uses a conditional Variational Auto-
Encoder model to generate grasps postures conditioned on additional variables
such as the grasp size. In this work we augment this model to also generate
grasp postures conditioned on the grasp type. The model is trained on a set of
labeled grasp samples acquired by teleoperating a robotic hand using a data-
glove. Using this model we are able to abstract away the low-level control of
each joint of each finger and generate grasps based on more general
characteristics such as the type and the size of the grasp. In this way we can
control all the fingers jointly by a single value, the grasp size, thus
greatly reducing the control parameters. In addition we are able to use the
same control algorithm for different precision grasp types, by changing the
grasp type conditional variable.

Finally, we can modify our controller to release objects instead of grasping
them. Given the pose of the hand in the world coordinate frame, which we can
acquire from the robotic arm that is attached to, we can use the forward
kinematics of the hand to compute the poses of each fingertip. Then using the
force readings of each fingertip we can calculate the global direction of the
net tangential force. If the angle between the direction of the net tangential
force and the direction of gravity is less than $90$ degrees, i.e. the net
tangential force’s direction is towards the ground, we assume that the
tangential force is due to gravity pulling the object, so the force controller
tries to grasp it. If the angle is more than $90$ degrees, i.e. the net
tangential force’s direction is upward, it means that something is pushing (or
pulling) the object upward, in which case we assume that the object is
touching on a support surface or someone is pulling the object so the
controller increases the grasp size given to the posture mapping function
proportionally to the normal force measured thus slowly releasing the object.
Opening the grasp is done by controlling the grasp size variable as follows:

$g_{size}(n+1)=g_{size}(n)+K\cdot f_{n}$ (2)

That way we can place objects on surfaces but also perform robot to human
handovers, where the robot holds the object and the human grasps the object
and slightly pushes or pulls it up, signaling to the robot that there is a
support surface. The robot then slowly releases the object by opening its
grasp. We showcase these scenarios in the experiments’ section.

Based on these observations, we present our force controller in Figure 2. The
hand starts in an open pre-grasp position, a latent point is sampled from the
prior distribution of the posture mapping function, and given the desired
grasp type and the grasp size a grasp posture, i.e. the joint angles of the
fingers, is sampled. The initial grasp size is set to the maximum value, and
when the force controller comes into effect and depending on the state of the
system and the forces on the fingertips grasp size changes by some value $C$,
according to equations 1,2, until the desired normal force is achieved.

To choose between grasping or releasing an object we use a finite state
machine formulation. When the hand reaches the desired grasp pose, which we
assume is provided, the GRASP state is activated, in which the controller
tries to grasp the object. When the controller detects that the tangential
force applied to the object is coming from a support surface the state changes
to the RELEASE state, in which the controller releases the object by opening
the grasp. You can see the full algorithm in Python-like pseudocode in Figure
3.

To summarize, the advantages of our controller compared with previous
approaches are threefold: 1) instead of controlling each joint of each finger
of the hand we use only two variables, the grasp size and the grasp type,
which allows us to perform multiple grasp types by changing only one variable
while the grasp size variable is common among all grasp types, that greatly
reduces the complexity of the control process compared to independently
controlling a 21 DoF hand to perform different grasp types, 2) we do not rely
on slip prediction for controlling the desired normal force, which involves
gathering labeled data and works only for the hand poses in the training
dataset, and 3) we can use our controller to also release objects instead of
only grasping them.

Figure 4:  Our first experiment. The robot picks up a bottle, transports it,
and places down on the desk. In the bottom part of the figure, you can see the
control signals during this task.

##  IV Experimental Results

###  IV-A Experimental Set-up.

Figure 5:  The household objects used in our experiments.

For our experiments we used the Seed Robotics RH8D Hand [17], which is a
robotic hand with 7 DoFs. The hand is equipped with the FTS-3 force sensors
[15] in each fingertip, which are high resolution tactile sensors that provide
the 3D force applied in each fingertip. The sensor provides data at a rate of
50Hz. For the experiments the hand was mounted on a Kinova Gen3 7DoF robot. To
train the posture mapping function we used the CyberGlove to teleoperate the
hand and collect 468 grasps belonging to three precision grasp types: tripod,
pinch, lateral tripod. The architecture of the cVAE model was the same as in
[2], with the addition of the grasp type as a conditional variable, which was
one-hot encoded. We used 10 household objects shown in Figure 5. With the
heaviest object weighing $380g$ and the lightest $1g$. During the experiments
the trajectories of the arm were prerecorded, while the hand was controlled
online by our control algorithm.

Figure 6:  In the upper row of images, you can see our second experiment. The
robot picks up the chips can, rotates it 90 degrees, and places back down. In
the middle row, for our third experiment, the robot picks up the chips can,
rotates it 90 degrees, and hands it over to a person. In the bottom row, for
our forth experiment, the robot picks up a foam brick, rotates it 180 degrees,
and hands it over to a person, using a pinch grasp. Figure 7:  In our fifth
experiment, a person hands over an empty plastic cup to the robot, throws
coins in it to increase its weight while the robot adjusts its grip to
stabilize the object, and then hand overs the cup back to the person.

###  IV-B Parameter tuning.

To select the values of the parameters in our controllers we conducted
preliminary experiments where we tested lifting and releasing several objects,
with different physical properties. To select the value of the normal offset
force $f_{n}^{offset}$, we used an empty plastic cup as our test object, and
we choose a value such that the fingers do not deform the cup. The final value
of the parameter was set to -50 mN. To select the values of the gain $G$ and
the rate of decrease $K$, of the grasp size, we experimented with the heaviest
object in our dataset, which is the mustard bottle and weighs $380g$. The gain
$G$ was set to $2.0$ such that the desired normal force would be enough to
hold the object. The rate of change of the grasp size was set to $100.0$,
based on the operating frequency of the force sensor and the range of values
of the tangential force. For the tangential force averaging process we used a
parameter value of $\alpha_{t}=0.7$, because we want the controller to be
sensitive to fast changes in its value, that can arise for example during
lifting an object. For the normal force averaging process we used a parameter
value of $\alpha_{n}=0.5$, as we do not want it to be affected by noise that
could make the controller overconfident.

###  IV-C Experiments.

To explore the capabilities of our controller, we demonstrate five experiments
of increasing complexity: 1) we picked and placed a bottle using a tripod
grasp, 2) we picked, rotated and placed a chips can on a box using a tripod
grasp, 3) we picked, rotated and handed over the chips can to a person using a
tripod grasp, 4) we picked, rotated and handed over a brown foam brick to a
person using a pinch grasp, 5) a person handed over a plastic cup to the
robot, filled it with coins to increase its weight, and the robot then handed
it back to the person using a tripod grasp.

You can see the execution of the first experiment in Figure 4. Under the
pictures of the execution you can see the signals recorded by the controller:
the average normal force applied by all fingers (blue line), the thresholds
$f_{n}^{threshold\\_high}.$ (purple dashed line) and
$f_{n}^{threshold\\_low}.$ (yellow dashed line), the average tangential force
(green), and the grasp size used in each time-step (red). The task is divided
four stages: 1) (red part) the initial grasp of the object, in this stage the
force controller closes the grasp until the applied normal force is below the
offset $f_{n}^{offset}$, 2) (green part) the robot lifts the object, as it
tries to lift the tangential force increases, increasing the threshold, so the
grasp size decreases to apply more normal force, 3) (orange part) the robot
transports the object, you can see, in point A in the Figure, a perturbation
in the tangential force when the robot begins to move, the controller responds
by decreasing the grasp thus stabilizing the object, and 4) (blue part) the
robot enters the releasing phase, where it lowers the arm until it detects
that the tangential force is due to a support surface, then it stops lowering
the arm and increases the grasp size slowly releasing the object. In point B
in the Figure, you can see that there is noise in the tangential force, due to
the arm moving to place the object on the table, that is also reflected in the
desired normal force. Because we use the desired normal force as a threshold
and not as a reference signal this noise is not manifested in the control of
the grasp size. You can see the execution of the second experiment in the
upper part of Figure 6. This experiment demonstrates the ability of the
controller to handle arbitrary hand poses. The experiment is divided in four
parts: 1) the robot enters the GRASP phase and the force controller generates
grasps to achieve a normal contact force below the $f_{n}^{offset}$ threshold,
2) the robot lifts the object and adjusts the grasp size to avoid the object
falling, 3) the hand rotates to place the chips can on the horizontal
position, and 4) the robot enters the RELEASE phase, and the arm lowers until
the object touches the box, when the hand detects the supporting surface, it
starts to slowly release the object. You can see the execution of the third
experiment in the middle part of Figure 6. This experiment demonstrates the
ability of the controller to perform robot to human handovers. The experiment
is divided in four parts: 1) the robot enters the GRASP phase and the force
controller generates grasps to achieve a normal contact force below the
$f_{n}^{offset}$ threshold, 2) the robot lifts the object and adjusts the
grasp size to avoid the object falling, 3) the hand rotates to place the chips
can on the vertical position, and 4) the robot enters the RELEASE phase, the
arm stays still, the human grasps the object from the bottom and slightly
pushes it up, the hand then detects that there is a supporting surface and
starts to slowly release the object. You can see the execution of the fourth
experiment in the bottom part of Figure 6. This experiment is similar to
previous one, but the grasp type that the robot uses is a pinch grasp, that
involves only the thumb and the index finger. To perform this we only had to
alter the grasp type conditional variable that was given to the posture
mapping function. You can see the execution of the fifth experiment in the
bottom part of Figure 7. In the first part (blue) of the experiment the robot
closes its grasp, by reducing the grasp size, until the normal force is below
the force offset. In the next three parts (pink, green, red) the person throws
coins in the cup to increase its weight. You can see in the signal plots that
each time coins are added the tangential force decreases so the normal force
threshold decreases too. The grasp sizes then decreases as well in order to
apply more normal force. This experiment demonstrates the ability of the
controller to handle perturbations in the weight of the object during
grasping.

##  V Conclusion

In summary, we presented a controller that uses force feedback integrated with
conditional synergies to control a dexterous robotic hand to grasp and release
objects. We demonstrated that our controller can lift objects of different
weights and materials while avoiding slip, react online when the weight of the
object changes, place them down on surfaces, and hand them over to humans. In
addition, the control architecture is modular, so the synergy grasp mapping
component can be easily changed in order to control several precision grasp
types. However, our experiments also revealed various limitations of our
controller. For example our method fails to stabilize the object when
rotational slip occurs. In addition hardware limitations such as, slow update
rates and noise in the force measurements can create problems that result in
the object falling. In future work we plan to incorporate additional sensing
modalities, such as vision to alleviate some of these issues.

## Acknowledgment

Work partially supported by the H2020 FET-Open project Reconstructing the
Past: Artificial Intelligence and Robotics Meet Cultural Heritage (RePAIR)
under EU grant agreement 964854, and by the Lisbon Ellis Unit (LUMLIS), and
the FCT PhD grant [PD/BD/09714/2020].

## References

  * [1] R. S. Johansson and J. R. Flanagan, “Coding and use of tactile signals from the fingertips in object manipulation tasks,” _Nature Reviews Neuroscience_ , vol. 10, pp. 345–359, 2009.
  * [2] D. Dimou, J. Santos-Victor, and P. Moreno, “Learning conditional postural synergies for dexterous hands: A generative approach based on variational auto-encoders and conditioned on object size and category,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 4710–4716.
  * [3] H. Yousef, M. Boukallel, and K. Althoefer, “Tactile sensing for dexterous in-hand manipulation in robotics-a review,” _Sensors and Actuators A-physical_ , vol. 167, pp. 171–187, 2011.
  * [4] G. Westling and R. S. Johansson, “Factors influencing the force control during precision grip,” _Experimental Brain Research_ , vol. 53, pp. 277–284, 2004\.
  * [5] J. R. Flanagan, M. C. Bowman, and R. S. Johansson, “Control strategies in object manipulation tasks,” _Current Opinion in Neurobiology_ , vol. 16, pp. 650–659, 2006.
  * [6] N. Wettels, A. R. Parnandi, J.-H. Moon, G. E. Loeb, and G. S. Sukhatme, “Grip control using biomimetic tactile sensing systems,” _IEEE/ASME Transactions on Mechatronics_ , vol. 14, no. 6, pp. 718–723, 2009.
  * [7] G. De Maria, P. Falco, C. Natale, and S. Pirozzi, “Integrated force/tactile sensing: The enabling technology for slipping detection and avoidance,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , 2015, pp. 3883–3889.
  * [8] R. Krug, A. J. Lilienthal, D. Kragic, and Y. Bekiroglu, “Analytic grasp success prediction with tactile feedback,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , 2016, pp. 165–171.
  * [9] H. N. Dang and P. K. Allen, “Stable grasping under pose uncertainty using tactile feedback,” _Autonomous Robots_ , vol. 36, pp. 309–330, 2014.
  * [10] Z. Su, K. Hausman, Y. Chebotar, A. Molchanov, G. E. Loeb, G. S. Sukhatme, and S. Schaal, “Force estimation and slip detection/classification for grip control using a biomimetic tactile sensor,” _2015 IEEE-RAS 15th International Conference on Humanoid Robots (Humanoids)_ , pp. 297–303, 2015.
  * [11] F. Veiga, B. B. Edin, and J. Peters, “In-hand object stabilization by independent finger control,” _ArXiv_ , vol. abs/1806.05031, 2018.
  * [12] Z. Deng, Y. Jonetzko, L. Zhang, and J. Zhang, “Grasping force control of multi-fingered robotic hands through tactile sensing for object stabilization,” _Sensors (Basel, Switzerland)_ , vol. 20, 2020.
  * [13] B. A. Kent and E. D. Engeberg, “Robotic hand acceleration feedback to synergistically prevent grasped object slip,” _IEEE Transactions on Robotics_ , vol. 33, pp. 492–499, 2017.
  * [14] B. Siciliano and O. Khatib, Eds., _Springer Handbook of Robotics_. Springer Berlin Heidelberg, 2008. [Online]. Available: https://doi.org/10.1007%2F978-3-540-30301-5
  * [15] “Seed Robotics - RH8D Adult Robot Hand,” https://www.seedrobotics.com/rh8d-adult-robot-hand.
  * [16] H.-H. Ko, J.-S. Kim, J. Kim, J.-G. Baek, and S.-S. Kim, “Intelligent adaptive process control using dynamic deadband for semiconductor manufacturing,” _Expert Systems with Applications_ , vol. 38, no. 6, pp. 6759–6767, 2011\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0957417410014363
  * [17] “Seed Robotics - FTS3 3D High Resolution Tactile (Pressure) sensor,” https://kb.seedrobotics.com/doku.php?id=fts:fts3_pressuresensor.


# A graphical, scalable and intuitive method
for the placement and the connection of biological cells

Nicolas P. Rougier111INRIA Bordeaux Sud-Ouest, Talence, France –
<EMAIL_ADDRESS>
###### Abstract

We introduce a graphical method originating from the computer graphics domain
that is used for the arbitrary and intuitive placement of cells over a two-
dimensional manifold. Using a bitmap image as input, where the color indicates
the identity of the different structures and the alpha channel indicates the
local cell density, this method guarantees a discrete distribution of cell
position respecting the local density function. This method scales to any
number of cells, allows to specify several different structures at once with
arbitrary shapes and provides a scalable and versatile alternative to the more
classical assumption of a uniform non-spatial distribution. Furthermore,
several connection schemes can be derived from the paired distances between
cells using either an automatic mapping or a user-defined local reference
frame, providing new computational properties for the underlying model. The
method is illustrated on a discrete homogeneous neural field, on the
distribution of cones and rods in the retina and on a coronal view of the
basal ganglia.

Keywords: Stippling, Voronoi, Riemann mapping, topology, topography, cells,
neurons, spatial computation, connectivity, neural networks.

######

  1. 1 Introduction
  2. 2 Methods
    1. 2.1 Distribution
    2. 2.2 Connection
    3. 2.3 Visualization
  3. 3 Results
    1. 3.1 Case 1: Retina cells
    2. 3.2 Case 2: Neural field
    3. 3.3 Case 3: Basal ganglia
  4. 4 Discussion

##  1 Introduction

The spatial localization of neurons in the brain plays a critical role since
their connectivity patterns largely depends on their type and their position
relatively to nearby neurons and regions (short-range or/and long-range
connections). Interestingly enough, if the neuroscience literature provides
many data about the spatial distribution of neurons in different areas and
species (e.g. [25] about the spatial distribution of neurons in the mouse
barrel cortex, [20] about the neuron spatial distribution and morphology in
the human cortex, [4] about the spatial distribution of neurons innervated by
chandelier cells), the computational literature exploiting such data is rather
scarce and the spatial localization is hardly taken into account in most
neural network models (be it computational, cognitive or machine learning
models). One reason may be the inherent difficulty in describing the precise
topography of a population such that most of the time, only the overall
topology is described in term of layers, structures or groups with their
associated connectivity patterns (one to one, one to all, receptive fields,
etc.). One can also argue that such precise localization is not necessary
because for some model, it is not relevant (machine learning) while for some
others, it may be subsumed into the notion of cell assemblies [13] that
represent the spatiotemporal structure of a group of neurons wired and acting
together. Considering cell assemblies as the basic computational unit, one can
consider there is actually few or no interaction between assemblies of the
same group and consequently, their spatial position is not relevant. However,
if cell assemblies allows to greatly simplify models, it also brings implicit
limitations whose some have been highlighted in [23]. To overcome such
limitations, we think the spatial localization of neurons is an important
criterion worth to be studied because it could induces original connectivity
schemes from which new computational properties can be derived as it is
illustrated on figure 2.

Figure 1: Stippling. According to Wikipedia333 Stippling Wikipedia entry at
https://en.wikipedia.org/wiki/Stippling, Stippling is the creation of a
pattern simulating varying degrees of solidity or shading by using small dots.
Such a pattern may occur in nature and these effects are frequently emulated
by artists. The pair of boots (left part) have been first converted into a
gray-level image and processed into a stippling figure (right part) using the
weighted Voronoi stippling technique by [29] and replicated in [27]. Image
from [27] (CC-BY license).

However, before studying the influence of the spatial localisation of neurons,
it is necessary to design first a method for the arbitrary placement of
neurons. This article introduces a graphical, scalable and intuitive method
for the placement of neurons (or any other type of cells actually) over a two-
dimensional manifold and provides as well the necessary information to connect
neurons together using either an automatic mapping or a user-defined function.
This graphical method is based on a stippling techniques originating from the
computer graphics domain for non-photorealistic rendering as illustrated on
figure 3.

Figure 2: Influence of spatial distribution on signal propagation. A. A
k-nearest neighbours (k=5) connectivity pattern shows mid-range connection
lengths in low local density areas (left part) and short-range connection
lengths in high density areas (right part). B. Shortest path from top to
bottom using a k-nearest neighbours connectivity pattern (k=5). The lower the
density, the shorter the path and the higher the density, the longer the path.
On the far left, the shortest path from top to bottom is only 6 connections
while this size triples on the far right to reach 19 connections. Said
differently, the left part is the fast pathway while the right part is the
slow pathway relatively to some input data that would feed the architecture
from the top. C. Due to the asymmetry of cells position, a signal entering on
the top side (materialized with small arrows) travels at different speeds and
will consequently reach the bottom side at different times. This represents a
spatialization of time. D. Due to the asymmetry of cells position, a signal
entering on the left side (materialized with small arrows) slows down while
traveling before reaching the right side. This represents a compression of
time and may serve as a short-term working memory.

##  2 Methods

Blue noise [31] is an even, isotropic yet unstructured distribution of points
[22] and has minimal low frequency components and no concentrated spikes in
the power spectrum energy [35]. Said differently, blue noise (in the spatial
domain) is a type of noise with intuitively good properties: points are evenly
spread without visible structure (see figure 3 for the comparison of a uniform
distribution and a blue noise distribution). This kind of noise has been
extensively studied in the computer graphic domain and image processing
because it can be used for object distribution, sampling, printing, half-
toning, etc. One specific type of spatial blue noise is the Poisson disc
distribution that is a 2D uniform point distribution in which all points are
separated from each other by a minimum radius (see right part of figure 3).
Several methods have been proposed for the generation of such noise, from the
best in quality (dart throwing [8]) to faster ones (rejection sampling [6]),
see [18] for a review. An interesting variant of the Poisson disk distribution
is a non isotropic distribution where local variations follow a given density
function as illustrated on figure 3 where the density function has been
specified using the image gray levels. On the stippling image on the right,
darker areas have a high concentration of dots (e.g. boots sole) while lighter
areas such as the background display a sparse number of dots. There exist
several techniques for computing such stippling density-driven pattern
(optimal transport [22], variational approach [7], least squares quantization
[19], etc.) but the one by [29] is probably the most straightforward and
simple and has been recently replicated in [27].

Figure 3: Centroidal Voronoi Tesselation. A. Voronoi diagram of a uniform
distribution (n=256) where black dots represent the uniform distribution and
white circles represent the centroids of each Voronoi cells. B. Centroidal
Voronoi diagram where the point distribution matches the centroid
distribution.

###  2.1 Distribution

The desired distribution is given through a bitmap RGBA image that provides
two types of information. The three color channels indicates the identity of a
cell (using a simple formula of the type $identity=256\times 256\times
R+256\times G+B$ for $0\leq R,G,B<256$) and the alpha channel indicates the
desired local density. This input bitmap has first to be resized (without
interpolation) such that the mean pixel area of a Voronoi cell is 500 pixels.
For example, if we want a final number of 1000 cells, the input image needs to
be resized such that it contains at least 500x1000 pixels. For computing the
weighted centroid, we apply the definition proposed in [29] over the discrete
representation of the domain and use a LLoyd relaxation scheme.

${\bf C}_{i}=\frac{\int_{A}{\bf x}\rho({\bf x})dA}{\int_{A}\rho({\bf x})}$

More precisely, each Voronoi cell is rasterized (as a set of pixels) and the
centroid is computed (using the optimization proposed by the author that allow
to avoid to compute the integrals over the whole set of pixels composing the
Voronoi cell). As noted by the author, the precision of the method is directly
related to the size of the Voronoi cell. Consequently, if the original density
image is too small relatively to the number of cells, there might be quality
issues. We use a fixed number of iterations ($n=50$) instead of using the
difference in the standard deviation of the area of the Voronoi regions as
proposed in the original paper. Last, we added a threshold parameter that
allows to perform a pre-processing of the density image: any pixel with an
alpha level above the threshold is set to the threshold value before
normalizing the alpha channel. Figure 4 shows the distribution of four
populations with respective size 1000, 2500, 5000 and 10000 cells, using the
same linear gradient as input. It is remarkable to see that the local density
is approximately independent of the total number of cells.

Figure 4: Non-uniform distribution (linear gradient). Different population
distribution (size of 1000, 2500, 5000 and 10000 cells) using the same linear
gradient as input have been computed. Each distribution has been split into
four equal areas and the respective proportion and number of cells present in
the area is indicated at the bottom of the area. The proportion of cells
present in each areas is approximately independent ($\pm$2.5%) of the overall
number of cells.

###  2.2 Connection

Most computational models need to define the connectivity between the
different populations that compose the model. This can be done by specifying
projections between a source population and a target population. Such
projections correspond to the axon of the source neuron making a synaptic
contact with the dendritic tree of the target neuron. In order to define the
overall model connectivity, one can specify each individual projection if the
model is small enough (a few neurons). However, for larger models (hundreds,
thousands or millions of neurons), this individual specification would be too
cumbersome and would hide any structure in the connectivity scheme. Instead,
one can use generic connectivity description [11] such as one-to-one, one-to-
all, convergent, divergent, receptive fields, convolutional, etc. For such
connectivity scheme to be enforced, it requires either a well structured
populations (e.g. grid) or a simple enclosing topology [26] such as a
rectangle or a disc. In the case of arbitrary shapes as shown on figure 5,
these methods cannot be used directly. However, we can use an indirect mapping
from a reference shape such as the unit disc and take advantage of the Riemann
mapping theorem that states (definition from [5]):

Riemann mapping theorem (from [5]). Let $\Omega$ be a (non empty) simply
connected region in the complex plane that is not the entire plane. Then, for
any $z_{0}\in\Omega$, there exists a bianalytic (i.e. biholomorphic) map $f$
from $\Omega$ to the unit disc such that $f(z0)=0$ and $f^{\prime}(z0)>0$.

Such mapping is conformal, that it, it preserves angles while isometric
mapping preserves lengths (developable surfaces) and equiareal mapping
preserves areas. [17] introduced a method to compute the Riemann mapping
function using the Szegö kernel that is numerically stable while [30]
introduced numerical methods for solving the more specific conformal Schwarz-
Christoffel transformation (conformal transformation of the upper half-plane
onto the interior of a simple polygon). Furthermore, a Matlab toolkit is
available in [12] as well as a Python translation
(https://github.com/AndrewWalker/cmtoolkit) that has been used to produce the
figure 5 that shows some examples of arbitrary shapes and the automatic
mapping of the polar and Cartesian domains.

Figure 5: Conformal mappings. Examples of conformal mappings on arbitrary
spline shapes using the conformal Riemann mapping via the Szegö kernel [17].
Top line shows conformal mapping of the polar domain, bottom line show
conformal mapping of the Cartesian domain.

However, even if automatic, this mapping can be perceived as not intuitive.
Provided the shape are not too distorted, we’ll see in the results section
that ad-hoc mapping can also be used.

###  2.3 Visualization

Having now a precise localization for each cell of each population, we have
several ways of visualizing the activity within the model. The most
straightforward way is to simply draw the activity of a cell at its position
using a disc of varying color (a.k.a. colormap) or varying size, correlated
with cell activity. This requires the total number of cells to be not too
large or the display would be cluttered. For a moderate number of cells, we
can take advantage of the dual Voronoi diagram of the cell position as
illustrated on figure 2, using a colormap to paint the Voronoi cell. Finally,
if the number of cells is really high, A two-dimensional histogram of the mean
activity (with a fixed number of bins) can be used as shown on figure 8C using
a bicubic interpolation filter.

##  3 Results

We’ll now illustrate the use of the proposed method on three different cases.

###  3.1 Case 1: Retina cells

The human retina counts two main types of photoreceptors, namely rods, and
cones (L-cones, M-cones and S-cones). They are distributed over the retinal
surface in an non uniform way, with a high concentration of cones (L-cones and
M-cones) in the foveal region while the rods are to be found mostly in the
peripheral region with a peak density at around 18-20∘ of foveal eccentricity.
Furthermore, the respective size of those cells is different, rods being much
smaller than cones. The distribution of rods and cells in the human retina has
been extensively studied in the literature and is described precisely in a
number of work [10, 1]. Our goal here is not to fit the precise distribution
of cones and rods but rather to give a generic procedure that can be
eventually used to fit those figures, for a specific region of the retina or
the whole retina. The main difficulty is the presence of two types of cells
having different sizes. Even though there exist blue-noise sampling procedures
taking different size into account [35], we’ll use instead the aforementioned
method using a two stage procedure as illustrated on figure 6.

Figure 6: Cones and rods distribution. A. The density map for cones placement
(n=25) is a circular and quadratic gradient with highest density in the
center. B. The density map for rods placement (n=2500) is built using the rods
distribution. Starting from a linear density, “holes” with different sized are
created at the location of each cone, preventing rods to spread over these
areas during the stippling procedure. C. Final distribution of cones and rods.
Cones are represented as white blobs (splines) while rods are represented as
Voronoi regions using random colors to better highlight the covered area.

A first radial density map is created for the placement of 25 cones and the
stippling procedure is applied for 15 steps to get the final positions of the
25 cones. A linear rod density map is created where discs of varying (random)
sizes of null density are created at the position of the cones. These discs
will prevent the rods to spread over these areas. Finally, the stippling
procedure is applied a second time over the built density map for 25
iterations. The final result can be seen on figure 6C where rods are tightly
packed on the left, loosely packed on the left and nicely circumvent the
cones.

###  3.2 Case 2: Neural field

Neural fields describe the dynamics of a large population of neurons by taking
the continuum limit in space, using coarse-grained properties of single
neurons to describe the activity [34, 33, 2, 9]. In this example, we consider
a neural field with activity $u$ that is governed by an equation of the type:

$\tau\frac{\partial u(x,t)}{\partial
t}=-u(x,t)+\int_{-\infty}^{+\infty}w(x,y)f(u(y,t))dy+I(x)+h$

The lateral connection kernel $w$ is a difference of Gaussian (DoG) with short
range excitation and long range inhibition and the input $I(x)$ is constant
and noisy. In order to solve the neural field equation, the spatial domain has
been discretized into $40\times 40$ cells and the temporal resolution has been
set to $10ms$. On figure 7A, one can see the characteristic Turing patterns
that have formed within the field. The number and size of clusters depends on
the lateral connection kernel. Figure 7B shows the discretized and homogeneous
version of the DNF where each cell has been assigned a position on the field,
the connection kernel function and the parameters being the same as in the
continuous version. The result of the simulation shown on figure 7B is the
histogram of cell activities using $40\times 40$ regular bins. One can see the
formation of the Turing patterns that are similar to the continuous version.
On figure 7C however, the position of the cells have been changed (using the
proposed stippling method) such that there is a torus of higher density. This
is the only difference with the previous model. While the output can still be
considered to be Turing patterns, one can see clearly that the activity
clusters are precisely localized onto the higher density regions. Said
differently, the functional properties of the field have been modified by a
mere change in the structure. This tends to suggest that the homogeneous
condition of neural fields (that is the standard hypothesis in most works
because it facilitates the mathematical study) is actually quite a strong
limitation that constrains the functional properties of the field.

Figure 7: Non-homogeneous discrete neural field. A. Turing patterns resulting
from a continuous and homogeneous neural field with constant and noisy input.
B. Turing patterns resulting from a discrete and homogeneous neural field with
constant and noisy input. White dots indicate the position of the cells. Mean
activity is compute from the histogram of cells activity using $40\times 40$
bins. C. Localized Turing patterns resulting from a discrete and non-
homogeneous neural field with constant and noisy input. White dots indicate
the position of the cells. Mean activity is computed from the histogram of
cells activity using $40\times 40$ bins.

###  3.3 Case 3: Basal ganglia

The basal ganglia is a group of sub-cortical nuclei (striatum, globus
pallidus, subthamalic nucleus, subtantia nigra) associated with several
functions such as motor control, action selection and decision making. There
exists a functional dissociation of the ventral and the dorsal part of the
striatum (caudate, putamen and nucleus accumbens) that is believed to play an
important role in decision making [24, 3, 21] since these two regions do not
receive input from the same structures. For a number of models, this
functional dissociation results in the dissociation of the striatum into two
distinct neural groups even though such anatomical dissociation does not exist
per se (see [14]). Without any proper topography of the striatal nucleus, it
is probably the most straightforward way to proceed. However, if each group
would possess its own topography, it would become possible to distinguish the
ventral from the dorsal part of the BG, as illustrated on figure 8 on a
coronal view of the BG. We do not pretend this simplified view is sufficient
to give account on all the intricate connections between the different nuclei
composing the basal ganglia, but it might nonetheless help to have better
understanding of the structure because it becomes possible to link external
input to specific part of this or that structure (eg. ventral or dorsal part
of the striatum). This could lead to differential processing in different part
of the striatum and may reconcile different theories regarding the role of the
ventral and the dorsal part.

Figure 8: Coronal view of the basal ganglia. A. Scalable Vector Graphic (SVG)
source file defining each structure in terms of border (solid black lines),
major and minor axis (dashed lines), input (red line) and output (blue line).
Local density is given by the alpha channel and structure identity is given by
the color. In this coronal view of the basal ganglia, the Caudate is red
(RGB=(0.83,0.15,0.15)), the GPe is blue (RGB=(0.12,0.46,0.70)) and the GPi is
green (RGB=(0.17,0.62,0.17)). B. Distribution of 2500 neurons respecting the
local density and structural organization (Caudate: 1345 cells, GPe: 884
cells, GPi: 271 cells). Neurons receiving input are drawn in red, neurons
sending output are drawn in blue. Each neuron possesses two set of
coordinates: one global Cartesian coordinate set and a local curvilinear
coordinates set defined as the distances to the major and the minor axis of
the structure the neuron belongs to. C. Mean activity histogram of the
different structures using 32x32 bins and a bi-cubic interpolation filter.
Each bin includes from zero to several neurons. D. Cell activities represented
using the dual Voronoi diagram of the cell position. Each Voronoi region is
painted according to the activity of the corresponding centroid (i.e. neuron).

##  4 Discussion

We’ve introduced a graphical, scalable and intuitive method for the placement
and the connection of biological cells and we illustrated its use on three
use-cases. We believe this method, even if simple and obvious, might be worth
to be considered in the design of a new class of model, in between symbolic
model and realistic model. Our intuition is that such topography may be an
important aspect that needs to be taken into account and studied in order for
the model to benefit from structural functionality. Furthermore, the proposed
specification of the architecture as an SVG file associated with the
scalability of the method could guarantee to some extent the scalability of
the properties of the model.

Notes: All figures were produced using the Python scientific stack, namely,
SciPy [16], Matplotlib [15] and NumPy [32]. All sources are available on
GitHub [28]

## References

  * [1] Peter K Ahnelt and Helga Kolb  “The mammalian photoreceptor mosaic-adaptive design”  In _Progress in Retinal and Eye Research_ 19.6 Elsevier BV, 2000, pp. 711–777  DOI: 10.1016/s1350-9462(00)00012-4
  * [2] Shun-Ichi Amari  “Dynamics of pattern formation in lateral-inhibition type neural fields”  In _Biological Cybernetics_ 27.2 Springer Nature, 1977, pp. 77–87  DOI: 10.1007/bf00337259
  * [3] B.. Balleine, M.. Delgado and O. Hikosaka  “The Role of the Dorsal Striatum in Reward and Decision-Making”  In _Journal of Neuroscience_ 27.31 Society for Neuroscience, 2007, pp. 8161–8165  DOI: 10.1523/jneurosci.1554-07.2007
  * [4] Lidia Blazquez-Llorca et al.  “Spatial distribution of neurons innervated by chandelier cells”  In _Brain Structure and Function_ 220.5 Springer Nature, 2014, pp. 2817–2834  DOI: 10.1007/s00429-014-0828-3
  * [5] Michael Bolt, Sarah Snoeyink and Ethan Van Andel  “Visual representation of the Riemann and Ahlfors maps via the Kerzman–Stein equation”  In _Involve, a Journal of Mathematics_ 3.4 Mathematical Sciences Publishers, 2010, pp. 405–420  DOI: 10.2140/involve.2010.3.405
  * [6] Robert Bridson  “Fast Poisson disk sampling in arbitrary dimensions”  In _ACM SIGGRAPH 2007 sketches on - SIGGRAPH ’07_ ACM Press, 2007  DOI: 10.1145/1278780.1278807
  * [7] Zhonggui Chen et al.  “Variational Blue Noise Sampling”  In _IEEE Transactions on Visualization and Computer Graphics_ 18.10 Institute of ElectricalElectronics Engineers (IEEE), 2012, pp. 1784–1796  DOI: 10.1109/tvcg.2012.94
  * [8] Robert L. Cook  “Stochastic sampling in computer graphics”  In _ACM Transactions on Graphics_ 5.1 Association for Computing Machinery (ACM), 1986, pp. 51–72  DOI: 10.1145/7529.8927
  * [9] “Neural Fields”  Springer Berlin Heidelberg, 2014  DOI: 10.1007/978-3-642-54593-1
  * [10] Christine A. Curcio, Kenneth R. Sloan, Robert E. Kalina and Anita E. Hendrickson  “Human photoreceptor topography”  In _The Journal of Comparative Neurology_ 292.4 Wiley-Blackwell, 1990, pp. 497–523  DOI: 10.1002/cne.902920402
  * [11] Mikael Djurfeldt, Andrew P. Davison and Jochen M. Eppler  “Efficient generation of connectivity in neuronal networks from simulator-independent descriptions”  In _Frontiers in Neuroinformatics_ 8 Frontiers Media SA, 2014  DOI: 10.3389/fninf.2014.00043
  * [12] Tobin A. Driscoll  “Algorithm 756, a MATLAB toolbox for Schwarz-Christoffel mapping”  In _ACM Transactions on Mathematical Software_ 22.2 Association for Computing Machinery (ACM), 1996, pp. 168–186  DOI: 10.1145/229473.229475
  * [13] Donald O. Hebb  “The Organization of Behavior: a Neuropsychological Theory”  New York: Wiley, 1949
  * [14] Mark D. Humphries and Tony J. Prescott  “The ventral basal ganglia, a selection mechanism at the crossroads of space, strategy, and reward.”  In _Progress in Neurobiology_ 90.4 Elsevier BV, 2010, pp. 385–417  DOI: 10.1016/j.pneurobio.2009.11.003
  * [15] J.. Hunter  “Matplotlib: A 2D graphics environment”  In _Computing In Science & Engineering_ 9.3 IEEE COMPUTER SOC, 2007, pp. 90–95  DOI: 10.1109/MCSE.2007.55
  * [16] Eric Jones, Travis Oliphant and Pearu Peterson  “SciPy: Open source scientific tools for Python”, 2001  URL: http://www.scipy.org
  * [17] Norberto Kerzman and Manfred R. Trummer  “Numerical conformal mapping via the Szegö kernel”  In _Journal of Computational and Applied Mathematics_ 14.1-2 Elsevier BV, 1986, pp. 111–123  DOI: 10.1016/0377-0427(86)90133-0
  * [18] Ares Lagae and Philip Dutré  “A Comparison of Methods for Generating Poisson Disk Distributions”  In _Computer Graphics Forum_ 27.1 Wiley-Blackwell, 2008, pp. 114–129  DOI: 10.1111/j.1467-8659.2007.01100.x
  * [19] S. Lloyd  “Least squares quantization in PCM”  In _IEEE Transactions on Information Theory_ 28.2 Institute of ElectricalElectronics Engineers (IEEE), 1982, pp. 129–137  DOI: 10.1109/tit.1982.1056489
  * [20] B.H. McCormick, R.W. DeVaul, W.R. Shankle and J.H. Fallon  “Modeling neuron spatial distribution and morphology in the developing human cerebral cortex”  In _Neurocomputing_ 32-33 Elsevier BV, 2000, pp. 897–904  DOI: 10.1016/s0925-2312(00)00258-7
  * [21] Matthijs AA Meer and A David Redish  “Ventral striatum: a critical look at models of learning and evaluation”  In _Current Opinion in Neurobiology_ 21.3 Elsevier BV, 2011, pp. 387–392  DOI: 10.1016/j.conb.2011.02.011
  * [22] Soham Uday Mehta, Brandon Wang and Ravi Ramamoorthi  “Axis-aligned filtering for interactive sampled soft shadows”  In _ACM Transactions on Graphics_ 31.6 Association for Computing Machinery (ACM), 2012, pp. 1  DOI: 10.1145/2366145.2366182
  * [23] Bhargav Teja Nallapu, Bapi Raju Surampudi and Nicolas P. Rougier  “The art of scaling up : A computational account on action selection in basal ganglia”  In _2017 International Joint Conference on Neural Networks (IJCNN)_ IEEE, 2017  DOI: 10.1109/ijcnn.2017.7965835
  * [24] J. O’Doherty  “Dissociable Roles of Ventral and Dorsal Striatum in Instrumental Conditioning”  In _Science_ 304.5669 American Association for the Advancement of Science (AAAS), 2004, pp. 452–454  DOI: 10.1126/science.1094285
  * [25] Joseph F. Pasternak and Thomas A. Woolsey  “The number, size and spatial distribution of neurons in Lamina IV of the mouse SmI neocortex”  In _The Journal of Comparative Neurology_ 160.3 Wiley-Blackwell, 1975, pp. 291–306  DOI: 10.1002/cne.901600303
  * [26] Hans Ekkehard Plesser and Håkon Enger  “NEST Topology user manual”, 2015  URL: http://www.nest-simulator.org/wp-content/uploads/2015/04/Topology_UserManual.pdf
  * [27] Nicolas P. Rougier  “[Re] Weighted Voronoi Stippling”  In _ReScience_ 3.1, 2017  DOI: 10.5281/zenodo.802285
  * [28] Nicolas P. Rougier  “Spatial Computation”, 2017  DOI: 10.5281/zenodo.1012306
  * [29] Adrian Secord  “Weighted Voronoi stippling”  In _Proceedings of the second international symposium on Non-photorealistic animation and rendering - NPAR ’02_ ACM Press, 2002  DOI: 10.1145/508530.508537
  * [30] Lloyd N. Trefethen  “Numerical Computation of the Schwarz–Christoffel Transformation”  In _SIAM Journal on Scientific and Statistical Computing_ 1.1 Society for Industrial & Applied Mathematics (SIAM), 1980, pp. 82–102  DOI: 10.1137/0901004
  * [31] Robert Ulichney  “Digital Halftoning”  Cambridge, MA, USA: MIT Press, 1987
  * [32] Stéfan Walt, S Chris Colbert and Gaël Varoquaux  “The NumPy Array: A Structure for Efficient Numerical Computation”  In _Computing in Science & Engineering_ 13.2 Institute of ElectricalElectronics Engineers (IEEE), 2011, pp. 22–30  DOI: 10.1109/mcse.2011.37
  * [33] Hugh R. Wilson and Jack D. Cowan  “A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue”  In _Kybernetik_ 13.2 Springer Nature, 1973, pp. 55–80  DOI: 10.1007/bf00288786
  * [34] Hugh R. Wilson and Jack D. Cowan  “Excitatory and Inhibitory Interactions in Localized Populations of Model Neurons”  In _Biophysical Journal_ 12.1 Elsevier BV, 1972, pp. 1–24  DOI: 10.1016/s0006-3495(72)86068-5
  * [35] Sen Zhang et al.  “Capacity constrained blue-noise sampling on surfaces”  In _Computers & Graphics_ 55 Elsevier BV, 2016, pp. 44–54  DOI: 10.1016/j.cag.2015.11.002