About the Camera Projection Model

This post is intended as a summary of my learning in the image processing using an RGB camera. My current workplace uses RGB cameras quite a lot in their engineering research. The brand I'm using is the Intel RealSense camera for the purpose of capturing human movement. The computation is done solely in Python using the "OpenCV" library.  Although image processing is a huge topic in itself, I'm interested merely in the calibration of the three-dimensional (3D) world system into the 2D camera image using OpenCV. Only one camera is assumed in this discussion.

1. Basic definitions
The way our human eyes see near things as bigger compared to those far away is called perspective projection. Whereas transformation is the transfer (or mapping!) of an object from one reference frame to another. When talking about space, we always talk about 3D space. On the other hand, an image produced by the camera is 2D. So overall, the perspective transformation deals with the conversion of the 3D world system into a 2D image.

In working with space, we need a so-called frame of reference. The 3D coordinate system has a 3D frame of reference. The basic principle of how human vision works is the same as how any camera works. Here we will use a pinhole camera concept, the simplest form of a camera. Refer to the figure. In the pinhole camera, the frame of reference commonly has an origin at the hole itself, which defines the center of projection of the camera coordinate system. In this way, the position of the pinhole camera in the camera coordinate system is at the origin. The 2D plane on which an image is clearly formed at the back of the camera is called the image plane. The plane is at a distance defined by the focal length from the projection center.

Fig-1: A simple representation of a pinhole camera.

2. The pinhole camera model

The pinhole camera model defines the perspective transformation from the camera coordinate system to an image coordinate system. 

The center of projection C is called the camera center or the optical center. The line connecting the camera center perpendicular to the image plane is called the principal axis or optical axis in Physics. The point where the principal axis passes through the image plane is called the principal point p. We assume that p is the origin of the image coordinate system on the image plane. An object can change its size on the image plane depending on whether it is near or far from the camera along the principal axis.

Fig-2: The pinhole camera model defines the mapping between the 3D coordinates into 2D coordinates. Look at the right panel. Suppose M is at (X, Y, Z) and m is at (x, y). Point C serves as the origin of the coordinate system. The proportionate triangle rule says that y / Y = f / Z; hence y = f .Y / Z. The same rule applies to x.



Let the center of projection C be the origin of the camera coordinate system, and the image plane is f away from point C. The image plane has been brought forward, noting that the position behind or in front of the camera center is congruent. Let M be a point in a 3D space in the camera coordinate system. We want to map this into a point m in the image plane. Using the proportionate triangle rule, the projection mapping from 3D space to 2D image coordinates is:

The above transformation above does not seem linear because of the division by Z. Let's introduce a new coordinate called the homogenous coordinate system. In this new coordinate, a new dimension is added in the coordinate such that any point (a, b) becomes (a, b, 1). Here the extra dimension of a point scales the existing dimensions. The above equation can then be expressed in the homogeneous coordinates format as a matrix multiplication:

3. The pixel coordinate system

 

Fig-3: Schematic shows forward projection from the real world coordinate to the pixel coordinate system.

The previous formula assumes that the origins of the camera and image coordinate systems are exactly along Z. In reality, this is not true (Fig-3). This is because an image is captured by a sensitive semiconductor device which has its own coordinate system. Although both the image and pixel coordinate systems are 2D, they are not exactly identical.

Suppose the principal point p has shifted, then the above formula can generalize into x = f.X/Z + px;  y = f.Y/Z + py.

Notice that the term focal length has changed. The f_x , f_y represent the focal length of the camera in terms of pixel dimensions. The pixel coordinate system is what the computer system hooked up with the camera will understand.

4. Representing the world system
In the equations previously, the subscript cam is written to denote that the points are represented in the camera's frame of reference. In general, a point in the 3D has its own frame of reference in the world coordinate system. For example, you can tilt either the camera in a certain direction, or the object, or both at the same time. We need another transformation matrix that maps the world in the camera reference frame. The transformation usually involves rotation and translation (like the affine transformation in the functional MRI data).

Recall that point M in Fig-2 is stated in the camera coordinate system. However, its position is not the position in the actual world coordinate system. Look at the figure below. The camera and world coordinate systems are related through a transformation matrix [ R | t ] = .

Fig-4: Representation of the transformation between the 3D world coordinate and 3D camera coordinate systems. The transformation matrix consists of R and t, where R denotes a 3×3 rotation matrix, and t is 3×1 translation matrix. The two matrices together defines the extrinsic parameters of a camera in the external world.



Accordingly, we can compute X_cam = R_x.U + t_x  ; Y_cam = R_y.V + t_y ; and  Z_cam = R_z.W + t_z .

Another way to look at the [ R | t ] matrix is the following. The R matrix represents 3 unit vectors (in each axis) of the camera coordinates in the world reference frame. The t, on the other hand, represents the origin of the camera coordinate system in the world reference frame.

Altogether, this leads us to the final transformation matrix. Suppose point M is located at (U, V, W) in the world coordinate system and you want to map this to a point m at a location (x, y) in the image coordinate system.

The above equation can be written as m = K [ R | t ] M. Here, let us define two more crucial terminologies:

1. Intrinsic parameter K, shows the internal orientation of a particular camera. It is fixed and unique for the camera, and is independent of the scene in the real world as long as the focal lengths do not change. 
2. Extrinsic parameters, R and t, show camera orientation or position to a world coordinate system. It translates coordinates of a point in the real world (world coordinate) to a coordinate system that is fixed with respect to the camera (camera reference).

This is the final equation for the forward projection. An inverse process, backward projection, retrieve the world coordinate system from the image position.

Camera calibration can be performed easily with the OpenCV package using a piece of checkerboard to obtain the transformation matrix. Often, this calibration also includes lens distortion correction inherent in any optical system.

Source:
1. HedVision Github
2. Calibration and 3D reconstruction

Looking at PCA in Action

There are different approaches to perform a PCA operation. Software tools (such as Matlab and R) offer two common approaches that we are going to discuss here.

1. Using Eigen-decomposition

The first approach is based on what was written in the previous blog post on PCA. The step-by-step example of PCA below is in Matlab as it is easier to work with matrices. The first step is to mean-center the original data. Suppose you have 2D matrix called data, with [x₁, x₂, ..., x₁₀₀] and [y₁, y₂, ..., y₁₀₀] with 100 observations.

% Step-1: Mean-center the original data.
data2(:,1) = data(:,1) - mean(data(:,1));
data2(:,2) = data(:,2) - mean(data(:,2));

Next, we find the covariance matrix to capture all possible relationships of the variables in the dataset.

% Step-2: Calculate covariance matrix of the mean-centered data.
C = cov(data2);

Once the covariance matrix C is computed we then perform eigen-decomposition on it to find the eigenvectors and eigenvalues according to CV = VD, leading to C = VDV^T. (I think we have to flip the matrix V left and right).

% Step-3: Produce eigenvalues and eigenvectors: a diagonal matrix D of eigenvalues and matrix V
% whose columns are the corresponding eigenvectors so that C*V = V*D.
[V,D] = eig(C);

Remember now the eigenvalues measure the "amount of variation" retained by each principal component. So, we rank the eigenvalues according to the magnitude or "importance". For example, suppose the columns of matrix V and the diagonal elements of D contain the following (see right):

The diagonal elements contain two eigenvalues λ₁ and λ₂. After ranking these two, it is seen that  λ₁ > λ₂, which means that the eigenvector v₁ corresponds to the first component, and the one that corresponds to the second is v₂.
Finally, the last step is to rotate or re-orient the data by [new data] = [original data].[eigenvector]^T. The new dataset is called the component scores.

% Step-4: The new axes is formed using V, to rotate the mean-centered data.
newdata = data2 * V';
newdata = fliplr(newdata) 

What is the percentage of information explained or accounted for by each component? It turns out that this can be found by dividing the corresponding eigenvalue by the total sum of all eigenvalues.

% Step-5: Get the variance accounted for, use matrix D containing the eigenvalues.
explained = D/sum(D(:))

Look at Fig-1. The original dataset in the left panel shows correlated variables. The middle panel shows how the dataset can be represented in terms of two principal components. PC1 axis is the first principal direction along which the data show the biggest variation. PC2 axis is the second most important direction and it is orthogonal to the PC1 axis. The new dataset in the right panel is a "projection" onto the principal components, where projection onto PC1 (X-axis) is uncorrelated with projection onto PC2 (Y-axis). 

Fig-1: The process of rotating the original data into a new dimension of principal components.

Where is the dimensionality reduction here? Take a look at the right panel above. The transformed data have just as many dimensions as the original data. By visual inspection, most variability is along the PC1 (X-axis) so you can just throw away the data plotted on the Y-axis. In other words, our 2D data can be reduced to 1D by projecting each sample onto the first principal component. Do we lose anything? Well, a little bit, depending on the variation that the PC1 and PC2 explain. In data science, data variation means information. Based on the above example, we find that PC1 and PC2 carry 1.284/(0.049+1.284) = 96%; and 0.049/(0.049+1.284) = 4% of the variance of the data respectively.

2. Using Singular Value Decomposition
PCA is intimately related to another statistical analysis called the Singular Value Decomposition (SVD). Although notations may differ in different textbooks, SVD essentially involves the factorization of matrix X = UDV^T, where U is an n × n matrix, with columns as orthogonal unit vectors of length n (left singular vector of X); D is an n × p rectangular diagonal matrix, the singular values of X; V is a p × p, with columns as orthogonal unit vectors of length p (right singular vectors of X).

The new data can be obtained through Y = XV = UD.  The transpose of V is sometimes called the whitening and can be used as a preparation in ICA. Columns of V multiplied by the square root of corresponding eigenvalues, i.e. eigenvectors scaled up by the variances, are the coefficients or loadings.

Now PCA can also be done in Matlab directly using a function:
[coeff,newdata,latent,tsquared,explained,mu] =  pca(data)

This Matlab function automatically mean-centers the data. The coefficients are also known as loadings. The loadings can be understood as the weights for each original variable when calculating the new scores. Each column of coeff contains coefficients for each component, and the columns are in descending order of the variance. If you type coeff and V in the Matlab terminal, you may find similarity. Recall that the eigenvectors are like the coefficients to predict the new scores from the original data (see Step-4 above!). But eigenvectors scaled up by the variances are called loadings. The terms can be confusing because the Matlab command already does normalization by default.

The newdata refers to the scores, which is a 2 × 100 matrix that represents the transformed data in the principal component space. The maximum number of principal components equals the number of variables in X. The component variances or latents are actually the eigenvalues of the covariance matrix of X. The latent values are arranged in descending order according to the magnitude of variance, i.e. ranking the eigenvalues. The function also returns explained, i.e. the proportion of the total variance explained by each component; and mu, the estimated mean of each variable in X.

To recover the original data, simply multiply the component scores with the transpose of matrix V. Remember to put back the mean value to each of the variables.

% Recover the original data
origdata = fliplr(newdata) * V' + mu;

In R, PCA can be performed automatically using a function princomp (using eigen-decomposition) or prcomp (using the SVD operation, for better accuracy in the results) from the built-in stats package. We will discuss one example in R below.

3. Case-study Using PCA
One nice package to easily extract and visualize the results of PCA is by using functions from the "factoextra" package that is based on "ggplot2". As the demo dataset, I am using the decathlon2 dataset from the "factoextra". Briefly, it contains a list of 27 athletes who have performed several different types of sports such as javelin throw, long jump, and so on. We call the athletes as individuals and the sports performance as the numeric variables. The dataset can be divided into:
- Active individuals (row. 1 to 23).
- Active variables (col. 1 to 10), to be used in PCA.
- Supplementary individuals (row. 24 to 27).
- Supplementary variables (col. 11 to 13), used in later prediction.

We first perform the PCA on the Active individuals/variables with the scaling option. This is an important step when you have values of different scales, e.g. metres and feet. After that, we can use the spree plot to visualize the different principal components ranked according to the proportion of variance explained. Recall that the eigenvalues reflect the amount of variance retained by the components. Eigenvalues are large for the first PC and small for the subsequent PCs. The first PC corresponds to a direction with the maximum amount of variation in the data.

We can ask other questions arising from the PCA results:

(1) How do we know the variance accounted after PCA?
Get this by extracting using the summary() function, we can see that the PC1 accounts for 41% of the total variance in the data, PC2 accounts for 18% variance, and so on. This is also what fviz_eig() provides. There are up to 10 original variables so we get up to 10 principal components.

(2) How many principal components are useful?
PCA does not tell us # components are useful, so there is no magic answer for this. From the cumulative proportion as shown above (or the spree plot), we can safely say that the first six components (PC1 - PC6) are indeed sufficient as they account for more than 90% of the total variance in the data.

(3) How are the original axes look like in terms of the components?
Before PCA, we have in the dataset the original axes which are formed by the 10 variables. After PCA, we get rotated axes which are formed by the principal components. To extract the rotated variables, we write get_pca_var(res.pca). Note that Dim.xx and PCxx are the same thing!

We can then plot how the original variables look like in terms of the first two components, PC1 and PC2. To visualize can be achieved using fviz_pca_var(), which will output the graph as shown below, Fig-2a on the left panel. The plot shows a correlation plot or relationships between all variables based on the first two PCs. Negatively correlated variables are positioned on opposite sides of the quadrants. The distance between variables and the origin measures the quality of the variables on the new map. Variables that are away from the origin are well represented on the map.

Fig-2: (a) Left: A correlation circle with a radius of 1.0 shows how the original variables relate to one another in the new dimensions of the principal components. (b) Right: a graphical representation of how each variable contributes to the principal component. It can also be seen that the representations of the original variables to PC7 to PC10 are minimal.

(4) What is the representation of each variable to the principal components?

We can use the function from the "corrplot" package to highlight the most contributing variables for each component. The function takes in var$cos2 as the argument. This cos2 shows the quality of representation of the variables on the component map (square cosine, squared coordinates).

Refer to Fig-2b. There are 2 ways to interpret this figure. First, it indicates a degree of representation of the original variables (the rows) on each principal component (the columns). Both panels in Fig-2 are somehow related. If a variable has a high cos2 value, it will be seen as a large and dark blue circle in Fig-2b, the same variable is positioned closer to the circumference of the plot in Fig-2a. A small circle indicates that the variable is not perfectly represented by the PCs. In this case, the variable is close to the center of the circle. Another way is to say how each principal component contributes to the variable. So, for a given variable, the sum of the cos2 on all the principal components is equal to one.


For further reading, see below:
(1) PCA in 6 steps
(2) Tutorial on PCA
(3) Making sense of eigenvector/eigenvalue

On Cluster-based Algorithms

I no longer post neuro/stroke-related blog, but more and more machine learning stuffs. Well, "Big Data" and AI/Machine Learning have entered a new beginning worldwide; it becomes a trend now. And this post is inspired by what I found from several online sources regarding clustering algorithms.

What and Why?
Generally speaking, cluster analysis means you group the dataset into different groups or clusters. It is commonly used to identify the underlying structures or patterns in the data. Cluster types are usually not known beforehand. Technically, any cluster-based algorithms try to maximize within-group similarity and maximize between-group distance. In machine learning, clustering is considered unsupervised learning, i.e. an algorithm that tries to discover unknown patterns in the dataset without any known reference (teacher), or prior knowledge of the labels. Lastly, clustering is related to a similar (but not the same!) algorithm called 'classification' which works on reference or labeled data through supervised learning processes. Look at the table below.

Cluster analyses have been widely used in data mining, biomedical research (genetic sequencing, phylogenetic, evolution), image processing, business analytics and social media, etc.

About k-means Algorithm
Let's talk about one of the most popular unsupervised learning algorithms out there, the k-means. This algorithm tries to partition the dataset into k distinct clusters (sub-groups), where each point or observation belongs to only one cluster. Before running the iteration, pre-processing ensures that the dataset belongs to the same scale on the coordinates. If these data points are not normalized, then it may lead to false grouping. Often, missing values are sometimes suboptimal to clustering.

Once your data is ready, we will start the k-means. As with any machine learning techniques, k-means involves an iterative process up to a termination point. In each iteration, it computes the cluster’s centroid, which is the arithmetic mean of all the data points that belong to that cluster.

Step1: Define the number of clusters, k.
Step2: Randomly select k different points to serve as the cluster centroid.
Step3: Assign all other data points to a cluster whose centroid is the nearest.
Step4: After that, calculate the new centroid position for each cluster.
Step5: Repeat Step 3 & 4 until the process reaches the termination criteria.

Objective: to minimize the sum of the squared within-cluster distance, that is, the
Euclidean distance between the data points and the cluster’s centroid.Termination criteria may include: centroids location do not change after some iterations, data points remain in the same cluster, or the user-defined max. number of iterations (e.g. 100) has been reached.

Refer to the GIF animation below. Each centroid is depicted as a white cross. Prior to grouping, data points are still black in color. After the grouping is done, the data points change color accordingly. In each iteration, k-means minimizes the total within-cluster distance which results in a new centroid location, and a shift in the color shade after. Indeed, some points can be reassigned to a different cluster in each iteration. The right panel shows what is known as the "Elbow method" (Fig-1).

Fig-1: Step-by-step process showing k-means clustering algorithm and its performance indicator (*).

In another example below, the initial guess produces three centroid locations in the left panel. The right panel shows the final outcome.

Fig-2: Initial and final centroid position of three different clusters (k = 3) after the iteration process stopped.

where w_ik= 1 for point x_i if it belongs to a certain cluster k; otherwise, w_ik= 0; and μ_k is the centroid of the x_i’s cluster. Minimizing this function contains two parts: first we differentiate J with respect to w_ik first and update the cluster assignments. Following this, we differentiate J with respect to μ_k, and recompute the centroids after the cluster assignments from the previous step.

How do we know whether our k-means has done a good job? Although a visual inspection is easy, one quantitative way is to see whether the # of clusters is already 'good enough' by using some evaluation metrics.

The first type is called the elbow-method, which uses the within-cluster sum of squared distance mentioned earlier. As shown in the figure above, the "elbow" here refers to the profound deflection downward that signifies a maximum improvement in the sum of squared distance. It is crucial to note that this metric continues to decrease monotonically as k increases.

The silhouette method, on the other hand, determines the degree of separation between clusters. It computes a certain coefficient whose value can take up anything in the interval [-1, 1]. Therefore, ideally the value should be as big as possible and closer to 1 to have a good clustering.

Hierarchical Clustering
In hierarchical clustering (HC), grouping is done in steps between two closest or most similar points, to produce certain hierarchy that is depicted as a dendrogram. There are essentially two techniques, the more popular agglomerative HC (bottom-up approach) and the divisive HC (top-down approach). Refer to the following Fig-3 for an illustration.

Fig-3: Two main types of hierarchical clustering, a different way of doing cluster-based analysis (tds com).

Refer to the nice GIF animation below with the different colors in each iteration. The vertical axis is a measure of similarity or closeness of either individual data points or clusters. The height of this vertical axis (Fig. 4) represents the Euclidean distance between the relevant data items grouped together.

Fig-4: Step-by-step process showing agglomerative HC algorithm and its dendrogram with colour-coding (*).

Suppose we work with the agglomerative HC:

Step1: Begin by treating each data point or observation as a separate cluster.
Step2: Compute a distance metric to identify the two points that are closest together.
Step3: Then, merge two points as a new cluster (linkage).
Step4: Repeat Step 2 & 3 for any pair of clusters obtained earlier until it reaches k clusters.

A set of distance or proximity metric between two data points/clusters is used, such as Euclidean distance, maximum distance, Manhattan distance, cosine distance, etc. This matrix is updated in each iteration to display the distance between the pair. 

The order in which clusters are joined is controlled by the linkage methods, meaning the manner you link the two clusters and categorize them hierarchically. There are generally 3 types of linkage in HC, the choice of which is up to us:
   (a) Single linkage: nearest-neighbour distance.
   (b) Complete linkage: furthest-neighbour distance.
   (c) Average linkage: taking the average distance.
   (d) Ward's method: taking the minimum sum of squared distance.

With the Ward's method, the sum of squares starts out at zero, because every point is in its own cluster. It then grows as more clusters are merged in the second iteration (each cluster now consists of more data). Given two pairs of clusters whose centers are equally far apart, Ward’s method will prefer to merge the smaller ones.
Refer to the diagram for more information on the agglomerative hierarchical clustering.


Final notes: Model-based approach
Now let's contrast the two algorithms. Clustering results are much more reproducible (consistent) in the HC if we were to repeat the analysis multiple times. Conversely, different results can be obtained in k-means since we started of with random centroids for every fresh analysis. Care should then be taken in choosing which algorithm and in interpreting the results. Another thing to note is that HC is computationally expensive and cannot handle big data well, but k-means clustering can. The time complexity of k-means is linear, i.e. O(n); while that of HC follows power rule, i.e. O(n²).

Although the k-means algorithm is popular, it suffers some drawbacks. First, the shape of the data is best to be spherical. It will be sub-optimal when the shape of data deviates from spherical shapes. Moreover, it will be confused when there are potentially two overlapping clusters as there is no obvious measure of uncertainty. This is because k-means is hard clustering, no way for a probabilistic partitioning. Lastly, although some methods exist to 'predict' the best number of clusters, it does not know the number of clusters from the data and requires it to be pre-defined.

Traditional clustering algorithms presented here are heuristic-based algorithms that derive clusters directly based on the data. In contrast, another form of clustering, model-based clustering attempts to address this concern and provide a soft assignment where observations have a probability of belonging to each cluster, hence, incorporating a measure of probability or uncertainty to the cluster assignments. For high-dimensional data, model-based approaches are preferred with some iterative methods called the Expectation-Maximization (EM). Unlike the k-means method which uses the Euclidean distance while calculating the distance between each point, the EM method uses more sophisticated statistical models, Gaussian Mixture Model (GMM), so making this a model-based approach. 

Briefly, the EM algorithm is often used to provide the functions more effectively. In the E-step, data points are assigned to the closest cluster according to the highest likelihood. In the M-step, the new centroid of each cluster is computed. In other words, assign the data point x_i to the closest cluster judged by its sum of squared distance from the cluster’s centroid. Iteration stops if the likelihood converges or stabilizes.

(*) Nice tutorial on k-Means
[*] Image source: Giphy website.

Stroke Rehabilitation (III): Impairments in Somatosensation

(1) Prevalence of somatosensory loss post-stroke

Evidence-based practice and research have established that impairments resulting from a stroke happen not only in motor domains (e.g. inability to perform reaching, loss of balance, and slurred speech) but also in somatosensory domains (loss of tactile sensation, limb position sense, and perceiving force). Carey LM (1995) stated that the loss of somatosensation occurs in about 60% of stroke survivors and has detrimental impacts on the quality of life, e.g. in spontaneous use of the hand and object manipulation. 

Intact somatosensation is essential for motor control since it is reliant on both intact feedforward and intact feedback from afferent inputs. It has been suggested that a learned non-use phenomenon with sensory loss leads to further deterioration of motor abilities. Despite this fact, the association between sensory impairment and outcomes following stroke has received limited focus in rehabilitation research. One reason is that most clinicians assume that spontaneous motor recovery occurs in the first 4 weeks post-stroke in the acute phase (Jia-Ching et al, 2005; Wing et al, 1990; Heller A et al, 1991; Lincoln NB et al, 1991), and that somatosensory recovery will arguably follow suit and therefore receives less attention. Another reason is the complexity in measuring the sensing ability of different modalities.

More recently, Connell et al (2008) conducted a newer prospective study with 70 patients with a first stroke assessed on admission day, and 2, 4, and 6 months post-stroke. Their findings did not contradict the earlier findings by Carey LM (1995). Of the sample collected, the authors found that 7–53% had impaired tactile sensations, 31–89% impaired stereognosis, and 34–64% impaired proprioception. Specifically, proprioception and stereognosis (the ability to perceive 3D shape and depth) were more frequently impaired than tactile sensations. This is in contrast to a study by Kim et al. [16] on acute stroke patients, who found that proprioception was less impaired when compared with localization and two-point discrimination regardless of the lesion.

Connell et al further said that the different somatosensory modalities showed only slight agreement between impairment within the same body areas, suggesting that the modalities are independent of each other. This suggests that it is necessary to include all somatosensory modalities while assessing one body part. In contrast, the high agreement between sensory modalities in adjacent body areas means that it is probably not necessary to assess all related body parts, e.g. there was redundancy between the wrist and hand, or between the ankle and foot.

(2) NSA and RASP scales

One challenge of sensory assessment post-stroke is the variety of sensing modalities of somatosensation, ranging from tactile or touch, pressure, position sense, movement direction, pain, to temperature. Another important barrier is the lack of standardization and low reliability of the clinical assessment scale (Winward, et al 1999). At the moment, there are three common clinical assessment scales for sensory impairment: the sensation parts of the Fugl-Meyer Assessment for UE/LE, Nottingham Sensory Assessment (NSA), and Rivermead Assessment of Somatosensory Performance (RASP). All sensory tests are conducted in the absence of vision.

NSA was developed as a standardized clinical sensory assessment, assessing both sides of the body and all areas. The original version uses a 5-scale rating system, has good intra-rater reliability (the same clinicians did multiple times), but poor inter-rater reliability (different clinicians did the same assessment) and was time-consuming (Lincoln NB, et al, 1991). The NSA measures tactile sensations (light touch, temperature, pinprick, pressure, tactile localization, and bilateral simultaneous touch), on the face, trunk, shoulder, elbow, wrist, hand, hip, knee, ankle, and foot, on both the paretic and normal side. The poor reliability has prompted a revision of NSA according to the Erasmus MC version (Em-NSA) (Stolk-Hornsveld, 2006) with lesser items to test but uses 3-scale rating system. Although NSA was shown to have concurrent validity with the more established and gold standard Fugl-Meyer Assessment for sensorimotor impairments (Scalha et al, 2011), this scale is still less attractive, with some clinicians view this to be a mere screening tool at best.

The RASP is a multi-modal sensory tool that tests six sensations (sharp/dull discrimination, surface pressure, tactile localization, temperature discrimination, joint movement, and joint movement direction discrimination), and two secondary sensations (sensory extinction and two-point discrimination) (Winward et al, 2002). The scale for proprioception was shown to have excellent test-retest reliability among sub-acute survivors and good concurrent validity with Motricity Index and Barthel Index, an ordinal scale used to measure performance in activities of daily living (ADL).

(3) Sensory impairments over time

One focus area in stroke rehab is the ability of a clinical assessment to predict recovery. The power of predictability helps clinicians to assess the stroke severity and to provide the most accurate intervention given a particular condition. Connell et al (2008) found the initial somatosensory impairment was significantly related to sensing ability at 6 months, accounting for 46–71% of the variance. The authors argued that the remaining factors were attributed e.g. to more cognitive factors, perceptual ability, and motivation. The spontaneous recovery over time was more obvious in the upper limb compared to the lower limb. 

In a study by Meyer and colleagues in Belgium (Meyer S. et al, 2016), the authors recruited acute stroke patients (< 1 week) and conducted sensory tests. Confirming earlier studies, they found that 41–63% of stroke survivors in the acute phase had a sensory loss in one of the modalities within the first week, but the deficits improved to 3–50% when assessed 6-month post-stroke. Proprioception score of Em-NSA moderately predicts motor ability at 6-month post-stroke as measured by the Fugl-Meyer UE and Action Research Arm Test. As a comparison, stereognosis moderately predicted motor ability at six months post-stroke as measured by Fugl-Meyer UE, the Motricity Index, and Action Research Arm Test. 

Using RASP, the somatosensory subtest of proprioception demonstrated the greatest level of recovery, but no patient achieved full recovery on all somatosensory subtests (Winward et al, 2007). In another study using the original NSA, Connell et al (2008) reported that most recovery of the upper limb tactile sensations and stereognosis occurred in the first 4 months, whereas recovery in proprioception continued over 6 months. Although the motor and functional recovery demonstrated continual improvement over time, somatosensory recovery showed marked variation in subtests both within and between patients. Both Connell's and Winward's studies reported that individual somatosensory modality recovery may be independent of other somatosensory modalities, despite the existence of parallel processing in the central nervous system.

Functional MRI has been beneficial to elucidate plastic changes in the brain following a stroke. Recently Carey et al. (2002) demonstrated in a single case study of severe sensory loss that there was little evidence of neural plastic changes in the early stages after stroke (2 weeks). However, a return of activation in ipsilesional primary and bilateral secondary somatosensory cortices was observed at 3 months and this was maintained at 6 months.

(4) Somatosensory recovery

Kessner et al (2016): Most cortical reorganization within the motor system occurs within 2–3 months post-stroke and stabilizes after 6 months. However, much less is known about the time course of recovery from somatosensory deficits, and about the mutual interaction of somatosensory and motor recovery. The deficits recovered at least partially, mostly within 3 months. Interestingly, some modalities (graphesthesia impairment and movement detection) even appeared to deteriorate during the time course after initial recovery (Julkunen, et al 2005). After about 6 months of recovery, the prevalence of somatosensory deficits seems to be lower. The initial somatosensory deficit was the strongest predictor for long-term somatosensory ability. Taken together, most patients recover at least partially from their somatosensory deficits, mainly during the first three to six months after stroke (Fig. 1). However, not all modalities necessarily recover with positive results. 

The probability and duration to reach a good rehabilitative plateau phase (Barthel index > 60, ambulation > 150 ft) are significantly worse, if stroke patients had combined motor and somatosensory deficits compared to motor deficits only. The authors summarized that somatosensory deficits after stroke have an important negative effect on motor and functional performance, especially proprioceptive impairments!

(5) Somatosensory-specific interventions

Retraining focusing on somatosensory-specific interventions is lacking. Stroke survivors with sensory impairments share that they often feel this is often a neglected aspect in their rehab program (Doyle, Bennet, Fasoli, & McKenna, 2010). 

A randomized controlled trial of sensory retraining has shown improvement in stroke patients’ sensory discrimination functions following a series of tactile and proprioceptive training (Carey et al., 2011). Such improvement was maintained and slightly increased at 6 weeks and 6 months follow-up. Evidently, the impact of somatosensory interventions on the recovery of sensation post-stroke is thought to be positive and significant. The authors have also argued that such a form of intervention is clinically beneficial, where reduction in deficit was targeted to improve lost abilities. 

On the other hand, Gopaul et al. (2018) proposed that the combination of somatosensory training with motor training may result in greater improvements in both motor and somatosensory functions, as compared to interventions that focus on individual function recovery. Most studies of such combined interventions have comprised active training components as it has the potential to drive neural plasticity (produce greater cortical activation extending to multiple areas), particularly when delivered insufficient dose. For example, Byl and colleagues (2008) demonstrated that patients who received higher-dose (72 hours) of integrated active-somatosensory and passive-motor training performed sensory discrimination tasks more accurately as compared to those receiving lower-dose (12-13.3 hours) of training. Despite the improvements observed, motor improvements were considered smaller than somatosensory improvements due to the reinforcement of passive movement training (mental practice and mirror therapy) in their study.

Basic Math for Principal Component Analysis

To begin, Principle Component Analysis (PCA) is a statistical tool for exploratory analysis that is famous for reducing the dimensionality of the data. It is a technique that involves vector space transformation to represent data in another form. It is useful for feature extraction and feature elimination, removing redundant features. The outcome of a PCA operation is a set of orthogonal variables, called the principal components. To learn PCA from a technical angle, I recommend the work by J Shlens (ref.[1]).

1. Intro to Matrix Operations
Before going into the math, there are a few technical terms in matrix algebra we have to be familiar with. First, eigenvector and eigenvalue. An eigenvector v is of a linear transformation (or simply a square matrix A) is a non-zero vector that changes only by a scalar factor (i.e. by the eigenvalue, λ) when that transformation is applied to it. We can say, [matrix].[eigenvector] = [eigenvalue].[eigenvector]; or we can write: Av = λv.  If the eigenvalue is negative, the direction is reversed. Refer to the figure on the right for a geometric description.

An eigenvector, with a real nonzero eigenvalue, points in a direction that is stretched by the transformation and the eigenvalue is the factor by which it is stretched.

Another term is about an orthogonal matrix. Let a square matrix A is invertible (it can be inversed) or non-singular.

A square matrix is orthogonal with real entries if, and only if, its transpose (A^T) is equal to its inverse (A^-1), that is AA^T yields an identity matrix I.  The square matrix is also normal if it satisfies A^TA = AA^T.

Importantly, if the matrix is orthogonal then it is necessarily invertible, unitary, and normal. The determinant of an orthogonal matrix is either +1 or −1. Next, we learn about another special matrix called a diagonal matrix, i.e. a square matrix with non-zero diagonal elements and all-zero remaining (non-diagonal elements). When a square matrix A is diagonalizable, there exists an invertible matrix B (a matrix that can have an inverse) such that B^-1AB produces a diagonal matrix.

Matrix S is symmetric only if it has a characteristic of S^T = S. Such matrix is also orthogonally (orthonormally) diagonalisable.

An example of a square matrix related to PCA is called the covariance matrix. More specifically, in this matrix we have the variances on the diagonal entries of this matrix, and the covariances on the off-diagonal entries. Take a look below. Here, the red box denotes the variance of variable height, whereas the green box shows the covariance between the score and age.

Suppose we have a huge dataset which has too much information. Can we re-express the original dataset optimally in a new form? Now the word 'optimality' means we do not want correlated or redundant features. After a series of matrix operations, the new or transformed dataset will have a new coordinate system. The whole idea of PCA is to force as much of the variation as possible into fewer dimensions, you can throw away the rest without losing much information. Let's relate PCA with the previous ideas on matrices:

In PCA, each eigenvector is a unit vector pointing in the direction of a new coordinate axis. The axis with the highest eigenvalue is said to be the axis that explains the largest variation, meaning an eigenvalue represents variance. The directions with the largest variances are the most “important” or principal. 

2. Technical Definition

Suppose we have a dataset of m × n matrix X, where the n columns are the number of observations and the m rows are the variables. We wish to linearly transform this matrix X into another matrix Y, also of dimension m × n, so that for some m × m transformation matrix P, we will obtain Y = PX.

So, PCA attempts to express correlated m variables into uncorrelated new variables through some sort of transformation. In matrix algebra, all possible correlations of any variable-pairs can be captured using the covariance matrix C_X. How do we decouple the variables? It turns out, we can get this by maximizing the variances or diagonal elements, and minimizing the co-variances or off-diagonal entries of the covariance matrix. So our goal is to obtain a new data Y such that the covariance matrix of the transformed data C_Y is strictly diagonal. To put it in another way, if we can find the transformation matrix P in such a way that C_Y is diagonal, then our objective is met.

The solution is based on an important property of matrix algebra on C_X called the eigen-decomposition. The effort to make C_Y a diagonal matrix is essentially to find P = E^T. How to get this new matrix? It turns out that the orthonormal matrix E can be found by diagonalizing the covariance matrix of the original dataset, C_X = EDE^T , where E is the matrix of eigenvectors which applies the transformation on X; and matrix D has diagonal elements containing the eigenvalues in descending order. In short, the principal components are formed by:
- The “new directions” of our data as indicated by the eigenvectors;
- The “magnitude”, or importance of each direction, as indicated by the eigenvalues.

Often, the normalized eigenvectors by the variance are also called the coefficients or loadings of the principal components. Remember, these coefficients define how the new principal component is formed as a linear combination of the original variables. Each column in the eigenvector matrix contains coefficients for each component in the descending order of the variance.

Dimensionality reduction works as follows. Suppose our original data is an m × n matrix X with m dimension and n data points.

If we want to transform the points to k dimensions (where, k < m), then select the first k eigenvectors of the matrix C_X sorted decreasingly according to the eigenvalues, and form a matrix with them, then use them as k × m matrix P. The resulting matrix is Y with size (k × m)(m × n) = k × n.

3. Closing remarks

PCA does not simply do a rotation. Normalization is always involved prior to PCA, which can be done by mean-centering the data and divide them with the standard deviation. This is important if our data contain variables of different scales, e.g. kilogram and pounds.

PCA and ICA are two slightly different things. While PCA expresses the data into a set of orthogonal variables, ICA simply represents data into an independent set of variables (not necessarily orthogonal). In ICA, the order of importance is irrelevant. Prior to ICA, data whitening is performed and PCA is one method to whiten the data. The figure on the right is common on the internet. PCA is predominantly used as a dimensionality reduction technique, for example, in computer vision, facial recognition, and image compression. It is also used for finding patterns in a high dimensional dataset in the field of finance, data mining, biomechanics, bioinformatics, etc.

References:
(1) "A tutorial on Principal Components Analysis" by J. Shlens (2014). Very nice technical tutorial!
(2) https://towardsdatascience.com/a-step-by-step-explanation-of-principal-component-analysis-b836fb9c97e2
(3) https://stats.stackexchange.com/questions/612/is-pca-followed-by-a-rotation-such-as-varimax-still-pca

Stroke Rehabilitation (II): Robotic Assessment of Proprioception

Previous studies have found the importance of somatosensory signals in motor control. Deaffarented patients are shown to have unique and distinct movement patterns compared to healthy control. Within the scope of motor behavior, movement-related somatosenses here means bodily senses that cover proprioception (sense of joint position), kinesthesia (sense of limb movement), and other senses originating from the mechanoreceptors, but does not include pain and thermal receptors. A large literature also shows that the somatosensory system is important for motor learning.

Neurological injuries (e.g. stroke, Parkinson's Disease, and multiple sclerosis) are often accompanied by somatosensory deficits. There has been a growing interest to assess and train the somatosensory system in a clinical setting. In the case of stroke, the somatosensory test is part of a standard clinical assessment. However, a recent systematic review (Connell & Tyson, 2012) has shown that such evaluation is known to be unreliable. The authors commented that the sensory section of the Fugl-Meyer Assessment and Erasmus version of Nottingham Sensory Assessment has the best balance between usability and robustness. In a stroke study, it was found that proprioception and stereognosis were more frequently impaired than tactile sensations and that there is little agreement in impairment among different somatosensory modalities (Connell et al, 2008). More research on assessments and interventions for somatosensory impairment post-stroke is warranted.

Broadly speaking, there are two classes in the assessment of proprioception: joint-position matching (JPM) tasks and psychophysical threshold method (PTM).

JPM tasks evaluate the ability to replicate the position or velocity of a joint angle or a limb position in the absence of vision. Participants respond by replicating the reference movement or position using the same limb (unilateral protocol) or the contralateral limb (bilateral). 

On the other hand, PTM tasks assess the ability to detect the intensity of proprioceptive stimuli, or discriminate two stimuli equal in nature but differ in amplitude/intensity. For example, the stimuli can be tactile, position, angle, or vibration. Participants respond typically with a verbal response.

Although some traditional tools such as goniometers can be used, the reliability of such measurements is unproven. For this reason, the use of robotic devices in the assessment of somatosensory integrity becomes popular. Robotics promise increased precision and accuracy, in addition to better reliability compared with standardized observer-based ordinal scales. Examples of such devices include KINARM (bimanual upper limb), MIT Manus (upper limb), Lokomat (lower limb), Wristbot (wrist), AMADEO (fingers). Other devices include an isokinetic dynamometer for the lower limb (Biodex).

Joint Position Matching Task
Joint position matching (JPM) tasks can be exercised on both lower and upper limbs, and either using unilateral or bilateral matching limbs. A thorough review of this test is described by Goble (2010). Instead of a pure sensory test, JPM is seen to involve some working memory component. Let's begin with the lower limb assessment. A few studies used Biodex lower limb devices to study lower limb proprioception. In a study by Willem et al (2002), subjects lied down in a supine position with the ankle in position 15° plantar-flexion. Active and passive joint-position sense using JPM method was assessed at the ankle, and muscle strength (isokinetic peak torque) was determined. This is also an example of a unilateral protocol, i.e. the same limb is used to replicate the movement.

Fig-1: An example of the joint position matching task using the same limb of lower extremities.

Another study investigates the hip and knee joint proprioception in healthy subjects and incomplete SCI subjects (Domingo and Lam, 2014) using a Lokomat exoskeleton. During the task, the subject's knee (or hip) was first moved passively by the robot into a position (joint angle θ). After a short break, the robot then moves it to a so-called distractor position. The subject was then asked to bring the joint back to the original position by using a joystick control (θ_S). The nominal value |θ ─ θ_S| is the matching error, reflecting the subject's performance. Refer to the figure above.

For the upper limb, studies mostly involved the more proximal limb, i.e. shoulder and elbow joints. A study with a bimanual protocol or two-arm matching task used two groups of subjects. In the first group, subjects moved both middle fingers to the same spatial location (extrinsic). The second group, however, was asked to move the right finger to a mirror-symmetric location of the left middle finger with respect to the body midline (intrinsic). The authors showed that bimanual accuracy is higher for tasks involving extrinsic coordinate (Iandolo, et al, 2015).

KINARM has proven to be a popular robotic tool for bimanual tasks studying sensorimotor behavior in the healthy and clinical population. Dukelow et al did a series of studies using KINARM in patients with stroke and see whether proprioception can be assessed objectively. In one study, 74 healthy subjects and 113 subjects with acute stroke (62 left-affected, 51 right-affected) performed a JPM task with vision occluded. The robot moved the most affected arm at a preset speed of 28 cm/sec, direction, and magnitude. Stroke patients were then asked to mirror-match the movement with their opposite, active arm (Semrau et al, 2013). For control subjects, both arms were tested. All subjects performed 36 total movements, for a total of 6 movements in each of 6 movement directions. See the figure below. The authors found that most stroke patients (69% of left-affected; median, 28.0°; 49% of right-affected subjects; median, 22.1°) made significantly larger directional errors than 95% of the controls (95% of range, <22.4°; median, 14.8°).

Fig-2: An example of the joint position matching task using bilateral limbs of upper extremities.

Recently, Konczak, Masia, and colleagues (2016) used a wrist robotic device and found the anisotropic nature of sensory acuity using JPM method on the dominant arm. Active matching acuity for flexion/extension was found to be 4.64 ± 0.24°; abduction/adduction: 3.68 ± 0.32°; supination/pronation: 5.15 ± 0.37°. A similar study was also performed with young children (Marini F, et al., 2017) and it was found that the anisotropic nature of sensory acuity does not change significantly across age. See the figure below for more illustration.

Fig-3: An example of the joint position matching task of the wrist (distal upper limb).

Psychophysical Threshold Method

Assessment using PTM yields a certain psychometric function with nominal value or threshold of detection in movement speed or position/angle. The method utilizes principles of Psychophysics, which is discussed in a separate blog entry on the "Human Sensory Perception". The shape of this function reflects the variability of responses about this threshold. The most common paradigm used in PTM tasks is a method of constant stimuli but one drawback of such method is the length procedure (Simo et al, 2014); and the same trial usually repeats until 3 to 5 correct judgments of the same stimulus are achieved. More recently, a revised and faster version of PTM paradigm was proposed by Mrotek et al (2017) using the "method of adjustment" during ten iterative trials. Here, subjects repeatedly adjust the magnitude of the stimulus until it is just perceptible, with an equal number of trials approaching that estimated threshold from below (i.e., starting from smaller stimulus magnitudes) and above (starting from larger magnitudes). Typical force magnitude and the threshold of detection can be found below for a hemiparetic patient (left) and normal control (right).

Fig-4: An illustration of a psychophysical method in assessing the threshold of detection.

Studies involving the lower limb have investigated the use of robotic devices in assessing knee joint proprioception. Using a custom-made device similar to the one by Biodex, Hurkmans et al. (2007) assessed the smallest detectable angular change in the knee joint. A similar method was used to find the smallest detectable passive knee movement in both sagittal as well as the frontal plane (Cammarata, et al., 2011). This method lets the subjects respond when they are just able to detect a change in joint position. In another study by Lam at UBC (Chrishold et al, 2016) in patients with spinal cord injury, Lokomat was used to test hip and joint again. Here, the joint was moved passively in 4 different movement speeds (0.5, 1.0, 2.0, and 4.0 deg/s) in both flexion and extension directions. The subject had to respond when the movement was first felt.

Assessment of proprioception in conjunction with motor adaptation has been studied by Ostry et al (2010). This study involves MIT-Manus robotic arm that produces a velocity-dependent curl field which perturbs the movement of the upper limb. Sensory acuity was assessed using an adaptive staircase method or PEST. The authors found that motor adaptation causes a shift in sensory acuity, a term called sensory recalibration (by another group, Henriques et al.). In another study, threshold detection of a wrist movement was assessed using a Wristbot by Konczak and colleagues. The authors found that the mean threshold for wrist flexion was 2.15°± 0.43° and 1.52°± 0.36° for abduction.

Conclusion
The robotic-based somatosensory assessment has gained popularity. It has great potential clinical applications in areas that involve human movements and movement disorders. Until now, there is no clinical biomarkers to predict somatosensory impairment in stroke (see: Boyd et al, 2018 for a review).


Reference: Based on a nice book chapter, "Robotic techniques for the assessment of proprioceptive deficits and for proprioceptive training" by Casadio et al., 2018. Figures are taken from the relevant individual citation therein.

Stroke Rehabilitation (I): Neurophysiology of Stroke Recovery

The topic of stroke rehabilitation is huge. In this first out of several blog posts, I will summarize the fundamental concepts of stroke. Subsequent posts will continue with more specific themes, e.g. intervention strategy, somatosensory relearning, and so on. 

Stroke can be classified into either hemorrhagic or ischemic stroke. Hemorrhagic stroke is caused by a rupture or internal bleeding in any of the cerebral arteries. On the other hand, ischemic stroke is due to blockage of blood supply to the brain, causing cell death due to lack of oxygen (infarct). This type of stroke is much more common to occur in either large arteries (atherosclerosis) or small penetrating arteries (lacunar infarct), or cardioembolism (the heart pumps blocking materials up to the brain). Neuroimaging has been useful as a diagnostic tool.

As most patients survive the initial injury, the next biggest challenge is the management of long-term impairment, limitation of daily activities or disability, and reduced participation or handicap. The main focus of the rehabilitation post-stroke is the recovery of the impaired movement (by physiotherapists), and the associated functions in daily living (by occupational therapists). Technically, we have to understand the difference between the recovery of function ("I am able to use the hand and arm in daily activities again") and resolution of impairment ("I can gain back my strength and movement"). Regardless of which, motor recovery after stroke is confusing. The term cannot be separated from compensatory mechanisms, where a new type of movement is produced to achieve the natural way prior to the stroke. Often during the assessment, clinicians do not separate motor compensation and recovery. In doing a prognosis, it has been known that there is a large inter-individual variability, patient heterogeneity. While clinical assessment has been very crucial in the early phase, there also exists a degree of inter-rater variability of the clinicians.

Stroke recovery in the early phase correlates to the resolution of dying tissue, edema, and inflammation (Furlan et al., 1996; Stinear and Byblow, 2014), while later recovery relates mainly to disinhibition of redundant neural circuits, recruitment of functionally homologous pathways, and the creation of neural connections to overtake the previous functions of the damaged neurons (Rossini et al., 2007; Murphy and Corbett, 2009; Ackerley et al., 2011). Interestingly, such processes may occur on the ipsilesional and contralesional hemispheres and are not completely understood (Hoyer and Celnik, 2011; Buetefisch, 2015).

Several studies have shown that patients commonly demonstrate increased M1 excitability on the contralesional hemisphere (equivalent to the ipsilateral hemisphere for healthy patients) for movements with the affected side (e.g. Shimizu et al., 2002; Butefisch et al., 2008; Murase et al., 2004; Ward and Cohen, 2004). Such theory, known as the interhemispheric competition model, says that an interhemispheric imbalance occurs in patients where the ipsilesional M1 no longer inhibits the contralesional hemisphere and the contralesional side appears to inhibit the ipsilesional, possibly through the transcallosal fibers. The magnitude of such an imbalance appears to positively correlate with the degree of motor impairment (Murase et al., 2004), and the interhemispheric imbalance in other functional networks may also contribute toward other cortical functional disruptions including neglect and aphasia.

Neuroimaging studies have shown that bilateral activation of the motor cortex leads to poorer motor recovery in most stroke patients. Conversely, studies by Nick Ward and colleagues show that a shift from bilateral activation to unilateral activation is a sign of good recovery. In other studies, it has been shown that this statement may have a limitation, i.e. significant mirror movements of the unaffected hand, causing an increase in contralesional activity.




In 2008, Krakauer's team studied the recovery of motor impairment using improvement in the Fugl-Meyer scale (FM) of the upper limb (UL). They defined recovery of impairment as the difference between FM score in the few days after stroke and at a later time point (3 months) (Prabakharan S. et al, 2008). They pointed to the idea of the proportional recovery rule which reflects spontaneous recovery. The maximum FM score for UL = 66. The proportional recovery rule states that, at 3 months, patients should get approximately 70% of their maximum potential recovery back. Example: a patient with moderate hemiparesis of 46 will recover (66-46) x 0.7 = 60. This rule has been validated in subsequent studies. Interestingly, some severe patients do not follow this rule while some other severe patients do. Such a categorical phenomenon bears two consequences. First, perhaps the current rehab therapy has little or, if it does, limited impact on the recovery within 3 months after the onset of stroke. Second, there are some underlying neurophysiological mechanisms unique to those non-fitter severe patients (see Krakauer et al, 2015).

TMS may provide a valuable assessment tool early in stroke. It can be used to test the functional integrity and excitability of the descending corticospinal pathways. Studies have shown that the ability to elicit MEP (motor evoked potential) within 2 weeks after stroke indicates a good corticospinal tract (CST) and it serves as a good predictor for recovery. However, most of these studies only targeted the upper limb, not the lower limb (LL), to the difficulty in accessing the "leg" area in the medial wall of the central sulcus. See Bembenek J.P, et al, (2012) for a review. Byblow et al (2015) did an important study using TMS very early in stroke, where they show that non-fitter patients do not follow the proportional recovery rule because they do not have intact CST useful for recovery. Patients who have their posterior limb of the internal capsule above 0.15, a threshold or point of no return, do not recover at later assessments. The resting motor threshold (RMT) also displays the proportional recovery rule. Interestingly, the authors show that adding regular therapy session does not yield significant results, suggesting that a more intensive behavioral intervention (e.g. using robots) may be needed. Feng W, et al (2015) did another relevant study where they formulated a neuroimaging biomarker of stroke. The authors used a weighted CST load, a method to better estimate the integrity of the CST in the ipsilesional hemisphere. An initial assessment using such measure, instead of an initial FM score, would be a more graded and sensitive predictor of recovery.