LISA: a MATLAB package for Longitudinal Image Sequence Analysis

TL;DR

LISA is a MATLAB package designed for comprehensive analysis of large longitudinal image sequences, offering automated and manual processing, registration, and statistical comparison tools to identify significant differences across images.

Contribution

The paper introduces LISA, a MATLAB toolkit that streamlines the processing, registration, and statistical analysis of longitudinal image sequences, addressing challenges like heterogeneity and noise.

Findings

Automated noise reduction using Gaussian mixture models.

Effective registration methods for contrasting images.

Generation of statistically significant difference maps.

Abstract

Large sequences of images (or movies) can now be obtained on an unprecedented scale, which poses fundamental challenges to the existing image analysis techniques. The challenges include heterogeneity, (automatic) alignment, multiple comparisons, potential artifacts, and hidden noises. This paper introduces our MATLAB package, Longitudinal Image Sequence Analysis (LISA), as a one-stop ensemble of image processing and analysis tool for comparing a general class of images from either different times, sessions, or subjects. Given two contrasting sequences of images, the image processing in LISA starts with selecting a region of interest in two representative images, followed by automatic or manual segmentation and registration. Automatic segmentation de-noises an image using a mixture of Gaussian distributions of the pixel intensity values, while manual segmentation applies a user-chosen intensity cut-off value to filter out noises. Automatic registration aligns the contrasting images based on a mid-line regression whereas manual registration lines up the images along a reference line formed by two user-selected points. The processed images are then rendered for simultaneous statistical comparisons to generate D, S, T, and P-maps. The D map represents a curated difference of contrasting images, the S map is the non-parametrically smoothed differences, the T map presents the variance-adjusted, smoothed differences, and the P-map provides multiplicity-controlled p-values. These maps reveal the regions with significant differences due to either longitudinal, subject-specific, or treatment changes. A user can skip the image processing step to dive directly into the statistical analysis step if the images have already been processed. Hence, LISA offers flexibility in applying other image pre-processing tools. LISA also has a parallel computing option for high definition images.

Figures

Tables

filename1	Name of the csv file that contains the first sequence of images
filename2	Name of the csv file that contains the second sequence of images
scantype	Data scanning method. Possible options include
	‘row’ : using a row identifier (required arguments : nrow, rowid)
	‘col’ : using a column identifier (required arguments : colid)
	‘no’ : using a blank row
nframe	Number of frames in each sequence
nrow	Number of rows in one frame or the height of the image
ncol	Number of columns in one frame or the width of the image
rowid	A characters string that identifies the start of a frame
colid	An integer number that identifies the column where the frame numbers are stored
preprocess	Logical value to inform if the images have been pre-processed. Possible options include
	‘T’ : True
	‘F’ : False (default)
segtype	Type of image segmentation. Possible options include
	‘auto’ : De-noises the image by choosing a cutoff value to divide noise to the region of interest by using a mixture of normal distributions of the intensity values (default)
	‘manual’ : De-noises the image through user-chosen cutoffs for the intensity values to divide noise to the region of interest
registtype	Type of image registration. Possible options include
	‘auto’ : Rotates and/or shifts images based on a mid-line found by a linear regression (no user interface is required) (default)
	‘manual’ : Rotates and/or shifts images based on a line specified by two points selected manually by a user (the first click sets a reference point and the second click draws the reference line)
roi	Select region of interest after image registration. Possible options include
	‘T’ : True
	‘F’ : False (default)
display	Output option. Possible options include
	‘basic’ : Only displays the P-Movie (movie of multiple comparison adjusted p-values) (default)
	‘all’ : Displays O-Movies (movies of the original two sequences of images), R-Movies (movies of the registered sequence of images), D-Movie (movie of the curated difference), S-Movie (movie of the non-parametrically smoothed differences), T-Movie (movie of variance-adjusted smoothed differences), and P-Movie
dimenpmap	Select the dimension of P-Movie. Possible options include
	2 : Two dimensions (default)
	3 : Three dimensions
parcomp	Parallel computing. Possible options include
	‘T’ : True
	‘F’ : False (default)

Equations

$$\tilde{Y}(i,j)=\{\begin{array}{ll}Y(i,j)& \text{if Y(i,j) > C}\\ 0& \text{otherwise}\end{array}.$$

$$f(y)=\sum _{g=1}^G\frac{{\alpha }_g}{{\sigma }_g}\varphi \left(\frac{y-{\mu }_g}{{\sigma }_g}\right)\equiv {\alpha }_1\cdot f_1(y)+(1-{\alpha }_1)\cdot f_2(y)$$

$$AvgCo{r}_j=\frac{1}{n-j}\sum _{i}cor(S_i^1,S_{i+j}^2)$$

$$\text{R}=\left[\begin{array}{ccc}cos\theta & -sin\theta & s_x\\ sin\theta & cos\theta & s_y\\ 0& 0& 1\end{array}\right]$$

$$d_j=m(x_{1j},x_{2j})+{ϵ}_j$$

$$\hat{\mathbf{\beta }}$$

$$\hat{m}(s)=\sum _{j=1}^N{p}_j(s)y_j$$

$$T(s)=\frac{\hat{m}(s)}{\hat{S}(s)}$$

Full text

K****eywords longitudinal images  $\cdot$ high-dimensional data  $\cdot$ image sequence  $\cdot$ multiple comparison  $\cdot$ non-parametric smoothing  $\cdot$ MATLAB

1 Introduction

Rapid development in imaging technologies has made medical imaging one of the most important techniques in providing essential evidence for clinical diagnosis and medical intervention. We can now collect unprecedented amounts of image data in both spatial and temporal dimensions per patient, from one or more modalities. The challenges in extracting valuable information from these images motivate modern updates and developments in image processing and analysis methods and tools. For example, three recent imaging books by Russ and Neal[1], Birkfellner[2], and Burger and Burge[3] show the fundamentals of image processing methods developed in the last score. Industries have also invested in developing image processing tools to accommodate many different modalities with types of images ranging from document imaging [4] and multimedia imaging [5] to medical imaging [6]. The National Institute of Health (NIH)’s Center for Information Technology is leading the Medical Image Processing, Analysis, and Visualization (MIPAV) project for quantitative analysis of medical images [7]. Most of the functionalities in MIPAV are for image (pre)processing of PET, MRI, CT, microscopy, 3D images and other widely-used images in medicine. One of the most widely used image analysis that MIPIV provides, uses summary statistics of the volume of interest (VOI) within a single image to generate contour VOI plots for all of the images to be analyzed.

Comparing images or assessing their similarities is one of the most important objectives in image analysis. The methods for such comparisons include the methods for pixel by pixel comparison and those based on a specific image measure [8]. The measure-based methods include Least Square Image Matching (LSM) [9], Hausdorff distance [10] and components of a matrix decomposition [11]. Among these measures, LSM has been used extensively because of its simplicity, but it may suffer from inaccurate initial input value and outliers [9]. To obtain comprehensive details of the differences over entire imaging region, a statistical pixel by pixel comparison is needed.

Denote an image sequence as spatial and temporal data $\boldsymbol{Y}(s,t,n)$, where $\boldsymbol{Y}$ is the intensity value at the pixel or spatial location $\boldsymbol{s}\in\boldsymbol{S}$ and time $\boldsymbol{t}\in\boldsymbol{T}$ for the subject $\boldsymbol{n}\in\boldsymbol{N}$. Here $\#\boldsymbol{S}$, $\#\boldsymbol{T}$, $\#\boldsymbol{N}$ are the sizes of sets $\boldsymbol{S}$, $\boldsymbol{T}$, and $\boldsymbol{N}$ respectively. Then, the number of measurements per subject, or the data dimension is $\boldsymbol{p}=\#\boldsymbol{S}\times\#\boldsymbol{T}$. This $\boldsymbol{p}$ can be huge, leading to a large-p data problem. If $\#\boldsymbol{N}$ is much smaller than $p$, we have a large-p-small-n problem. Large-p data challenges statistical multiple comparisons, no matter the number of subjects $\#\boldsymbol{N}$ is large or small. How to obtain an objective threshold for multiple comparisons of large, noisy and potentially heterogeneous images has challenged statisticians, especially for analyzing overflowing simultaneous images. Statistical parametric mapping (SPM) [12] is an effective and widely-used pixel by pixel image comparison approach. The motivation for developing Longitudinal Image Sequence Analysis (LISA) MATLAB package comes from generalizing and making a one-stop ensemble of image processing and analysis tool for complex, longitudinal images, by transforming the components of the Longitudinal Analysis and Self-Registration (LASR) procedure [13]. The LASR was developed to overcome the heterogeneity, and the large-p-small-n issue in assessing the effect of Neuromuscular Electrical Stimulation (NMES) experiments [14] in reducing the risk of developing pressure ulcers. The data for LASR were movies composed from sequences of seated pressure mat images obtained under a static and a dynamic stimulation, respectively, in each session from each visit by a patient, for eight patients observed up to 5 years. There were artifacts, heterogeneity, noises, and substantial spatial and temporal alignment issues from different sessions even on the same patient. The data was massive with a huge $p$ and a small $n$ of 8. These data required extensive imaging preprocessing and general statistical analysis method, Statistical Non-parametric Mapping (SNM) that is more general than the SPM, for each set of comparisons. Various code for the components of LASR was developed individually and specifically to the pressure data formate. It would have been easier if the LISA, a general purpose MATLAB package was available then. With increasing amounts of imaging data from different modalities in nowadays, it becomes imperative to implement and generalize the LASR method into a stand-alone MATLAB package for general application of images in standard image formate. Making LISA available to public as a one-stop MATLAB package sets a basis for future development for even large data and more generalized case. See the discussion in Section 6.

Figure 1 above shows the schematic structure of LISA tool, where the ‘image processing’ is called ‘image pre-processing,’ and ‘statistical analysis’ is called ‘image post-processing.’ The top part of Figure 1 presents the overall flow of LISA. A complete process of LISA goes through 2 steps: “I. Image Pre-processing” step as introduced in the left bottom flowchart of Figure 1 and “II. Image Post-processing (Statistical Analysis)” step as introduced in right bottom flowchart of Figure 1. However, if image pre-processing is unnecessary or a user chooses to use other MATLAB built-in tools to de-noise, segment or register their two sequences of images, LISA has the option to skip the image pre-processing step to go directly to image post-processing (statistical analysis). Image pre-processing is done by 1) selecting a region of interest (ROI), 2) de-noising the images by automatic or manual image segmentation, and 3) aligning the images by also automatic or manual image registration. Image post-processing is done automatically by 1) calculating curated differences, 2) de-noising the differences obtained from Step 1 by a non-parametric smoothing, 3) deriving variance-adjusted, smoothed values to get T-type statistics, and 4) providing final significances of the differences based on the multiplicity adjusted p-values. These four post-processing steps provide Statistical Nonparametric Mapping (SNM) with minimal assumptions on the images. Here we use the False Discovery Rate (FDR) to adjust the multiplicity.

The development of MATLAB package LISA(Longitudinal Image Sequence Analysis) has implemented and expanded the functionality of LASR to general cases for any images that can be aligned using a reference line (estimable from the image data directly) or general image processing tools in Matlab. A graphical user interface in LISA provides convenience for a user to quickly remove noise and artifacts, decide on the region of interests and choose a quick reference line if feasible. In addition, LISA provides a parallel computing option for a large sequence or high-resolution images.

In the remaining paper, we provide the LISA installation guide (§2), options for pre-processing (§3) and post-processing (statistical analysis) with a brief introduction of methodologies (§4), demonstrate the software LISA via an application (§5), and then conclude the article with a discussion, including limitations and further extensions (§6). For the details and interpretations of of the D, S, T and P maps, see Section 4. With features introduced in LISA, we improve the accuracy of statistical comparison for both static and dynamic sequences of images. LISA can tackle two sequences of time-series images in general, including seated pressure mat images used in the NMES study, time-series satellite image[15], and time-series chemical imaging data [16].

2 Installation

A step-by-step installation guide is provided in this section.

1. Download LISA App

LISA program “LISA.zip” can be downloaded from the following two locations.

$\cdot$

MathWorks File Exchange website : http://www.mathworks.com/matlabcentral/fileexchange/?term=type%3AApp

$\cdot$

Center for Statistical Research, Computing, and Collaboration(SR2c) website : http://sr2c.case.edu/software.html

The archive “LISA.zip” contains three files : LISAxxx.mlappinstall, LISAxxx.prj, and ReadMe where xxx is the package version number. “LISAxxx.mlappinstall” file is the installation file, “LISAxxx.prj"’ is a project file containing package information and “ReadMe” file is a quick starting guide.

2. Install and load LISA App

With a newer version of MATLAB (version 8 or above), one can install and load LISA by simply running the LISAxxx.mlappinstall file [17] and follow the three steps below. However, with an older version of MATLAB (version 7 or below) one can manually install LISA using the instructions provided in the following MathWorks website https://www.mathworks.com/help/matlab/matlab_env/manage-your-add-ons.html.

Run LISAxxx.mlappinstall file. Then MATLAB will be launched and a “Install” window will pop-up as shown in Figure 2 (a) 2. 2.

Continue by clicking the “Install” button. MATLAB will automatically install the package and will notify the installer from the “APP” tab when the installation is complete, as shown in Figure 2 (b) 3. 3.

Load LISA by clicking on the LISA App icon under “APPS - MY APPS” as shown in Figure 3

When LISA is loaded, a user can start using the program right away by running the lisa function. Typing “helpLISA” into MATLAB command window will provide a vignette including the syntax and examples of the lisa function.

3 Data and Image Processing

LISA is bulit to compare two sequences of images by contrasting each sequence with representing images taken in one session. The first stable image of each sequence is treated as its reference image, which will be used throughout the image processing. Data and image processing consists of the following 5 steps.

Data Scan : Scan image frames from “csv" files 2. 2.

Region of Interest(ROI) Selection : Crop images to the region of interest 3. 3.

Image Segmentation : Detect noise 4. 4.

Image Registration : Calibrate Images 5. 5.

Finer ROI Selection(optional)

3.1 Data Scanning

The current version of LISA accepts “csv" files as inputs. A user can convert various image formats to a “csv” format in MATLAB by reading images using imread and saving the image to a “csv” file using csvwrite. In general, the output of an imaging tool in a “csv” file includes information on the date, time, frame numbers, intensity values of the images and metadata such as equipment model, user information etc. Therefore, it is very important to distinguish the image intensity matrix in a numerical form to conduct statistical comparison. LISA is able to extract image intensity matrices from the raw data files that are in the following three formats as shown in figures 4(a), 4(b), and 4(c).

The data structure in Figure 4(a) requires a user to extract the data matrix using a row identifier. A user needs to set the parameter ‘scantype’ to ‘row’ and a row identifier character string needs to be passed over to the parameter ‘rowid’. LISA skips the header rows until it encounters the first identifier. It then scans next ‘nrow’ number of rows to extract as a numeric matrix and define it as the first image frame. LISA continues the same steps until it reaches the last identifier. Here is an example of the syntax for scanning the data using the row identifier having all other arguments with their default values.

lisa(‘file1.csv’,‘file2.csv’,‘row’,#frame,#row,#col,‘rowidentifier’,

‘’,‘’,‘’,‘’,‘’,‘’,‘’,‘’ )

The data structure in Figure 4(b) requires a user to extract the data matrix using a column identifier. A user needs to specify the column that contains information on the frame numbers by setting the parameter ‘scantype’ to ‘col’ and a column number needs to be passed over to the parameter ‘colid’. LISA reads the number in that column to extract numeric matrix and save each matrix to its corresponding frames. Here is an example of the syntax for scanning the data using the column identifier having all other arguments with their default values.

lisa(‘file1.csv’,‘file2.csv’,‘col’,#frame,#row,#col,‘ ’,# of column

identifier,‘’,‘’,‘’,‘’,‘’,‘’,‘’,‘’ )

The data structure in Figure 4(c) requires a user to extract the data matrix using a blank row right before the numeric intensity matrix. A user needs to set the parameter ‘scantype’ to ‘blank’ and LISA skips the header rows until it encounters the first blank row. It then scans next ‘nrow’ number of rows to extract as a numeric matrix and define it as the first image frame. LISA continues the same steps until it reaches the last blank row. Here is an example of the syntax for scanning the data using a blank row having all other arguments with their default values.

lisa(‘file1.csv’,‘file2.csv’,‘blank’,#frame,#row,#col,‘’,‘’,‘’,‘’,

‘’,‘’,‘’,‘’,‘’ )

3.2 Select Region of Interest

This is the first interactive step for a user. A user is given the opportunity to select a rectangular region of interest from the first images of each sequence. Cropping images to a specific area of interest enhances the result of data segmentation and data registration because this naturally eliminates some part of the noise. LISA uses the “imcrop" function within the MATLAB Image Processing Toolbox to select the region of interest.

3.3 Image Segmentation

It is critical to reduce noises and outliers by segmenting an image, which sets the noisy intensity values to zero before performing an image registration. LISA provides data-driven automatic and manual segmentation. A pixel is identified as noise if its intensity value is less than a threshold $C$. Let $Y(i,j)$ denote the pixel intensity value at the $i$th row and $j$th column of a representative image. Then the intensity values in the segmented image are $\widetilde{Y}(i,j)$ below:

[TABLE]

The threshold value $C$ is derived differently for automatic segmentation and manual segmentation, as shown below.

3.3.1 Automatic Segmentation

Automatic segmentation finds an optimal $C$ that separates data into two groups with the maximum separation, based on an estimated mixture of Gaussian distributions of the pixel intensity values, as introduced in [13]. The mixture density $f(y)$ has the following form:

[TABLE]

where $\phi$ is the standard normal density, $G$ is the number of mixtures, and $\{\alpha_{g}>0,\;g=1,...,G\}$ are mixing parameters satisfying $\sum_{g=1}^{G}\alpha_{g}=1$. The first-component density $f_{1}(y)=\phi\left[(y-\mu_{1})/\sigma_{1}\right]/\sigma_{1}$ represents the background distribution, while the second-component density $f_{2}(y)=\sum_{g=2}^{G}(\alpha_{g}/\sigma_{g})\cdot\phi\left[(y-\mu_{g})/\sigma_{g}\right]/(1-\alpha_{1})$ represents the distribution of pixel intensity inside the region of interest, also as a mixture of Gaussian distributions. The $C$ value automatically computed by LISA separates $f_{1}$ and $f_{2}$ at the valley bottom between the two densities - see the black vertical line in Figure 5.

Figure 5 is an example of the histogram and density plot for the pixel intensity values in a seated pressure map. The distribution of the intensity values is modeled by a mixture of 3-component Gaussian distributions. The threshold $C$ is determined by the vertical line that separates the first and second components in the fitted parametric Gaussian mixture model.

3.3.2 Manual Segmentation

Manual segmentation allows a user to choose $C$ visually based on a histogram of pixel intensity values. Figure 6 is a histogram of intensity values. It is natural to select the threshold $C$ at the gap where the target cross is located because at that threshold value the intensity values are distinctively divided into two groups. Once the threshold $C$ is selected, intensity values below $C$ will be replaced with zero, defined as the background noise.

3.4 Image Registration

There are two steps in image registration: temporal registration and spatial registration. Temporal registration excludes unstable images from each sequence and identifies the most similar subset of images to compare between the two sequences. Spatial image registration is done using a geometric transformation of the images.

3.4.1 Temporal Registration

The correlation coefficient of the intensities between two images is a reasonable measure of similarity between the two images. LISA temporally aligns two image sequences (of same speed) by shifting one sequence forward $j$ frames to match temporally the conditions of another sequence, using an Intensity-based Correlation Registration (ICR) below.

Exclude $5$ or the first 5$\%$ of images from each sequence, whichever the number is larger. This step removes noise images at the start of recording. Let $S^{1}_{1},...,S^{1}_{n}$ and $S^{2}_{1},...,S^{2}_{n}$ be the emaining $n$ sequential images in two image sequences $S^{1}$ and $S^{2}$. 2. 2.

Compute the shift-$j$ correlation coefficient of two sequences by averaging correlation coefficients of the data from the respective $i$-th and $(i+j)$th frames of two sequences over $i=1,...,n-j$, for each $j=0,...,n-1$ to obtain:

[TABLE]

Then find $j_{max}$ that maximizes the correlations coefficients $vgCor_{j}$ over $j$. In other words $j_{max}=\min\{j_{0}:AvgCor_{j_{0}}=\max_{j}(AvgCor_{j})\}$. 3. 3.

Align $S^{1}_{i}$ and $S^{2}_{i+j_{max}}$ for frames $i=1,...,n-j_{max}$.

3.4.2 Spatial Registration

LISA provides automatic registration and manual registration for the sequences after a temporal registration,

Automatic Registration

The following midline registration algorithm from Algorithm 3.1 in Wang et al.(2006)[13] is applicable to registering images that can be aligned by a rigid transformation using a midline (and a point) as a natural landmark.

Automatic Registration Algorithm:

Determine, column by column, the mid-points($m_{c}$) by counting non-zero values in each column of an image using $m_{c}$ = $\frac{r}{2}+\frac{(u-l)}{2}$, where $r,u,l$ are the number of rows, number of non-zero intensity values from the upper half and lower half of the column respectively for $c=1,...,\#ofcolumns$. 2. 2.

Fit a regression line through the mid-points and call this line the mid-line. 3. 3.

Conduct a rigid transformation by aligning each mid-line with a horizontal line as shown in Figure 7 using a rotation and location shift.

In other words, we define the transformation matrix R using the slope and intercept estimates of the mid-line from step 2 as

[TABLE]

where $tan\theta$ is the slope of the mid-line and $s_{x}$ and $s_{y}$ are the shifting parameters. Using $R$ we change $X$ to $X^{\prime}=RX$ where $X$ contains the information of the coordinate of the pixel defined as $X=[x,y,1]$.

Figure 8 shows an example of automatic registration. The image is translated to be aligned with the horizontal center line. The mid-line registration algorithm given above is fast to implement automatically when comparing sequences of many images. For images that require a rescaling or other types of registration, we suggest a user to use other MATLAB built-in image processing tools or manual registration in the next section.

Manual Registration

Manual registration aligns the images using a line defined by two user-selected points. Users select two points on the first images from each sequence which makes four point selections in total. The first points selected for each sequence are used as a reference point which is aligned at (5,# of rows/2) to shift the image to the same starting point. The second points selected is used to form a line which is used to rotate the image to the horizontal middle line of the image. An illustration of the point selection is provided in Figure 9.

Then LISA will perform a rigid transformation based on a shift and rotation matrix of R defined by the dotted blue line and 1st points from Figure 9, where R = $\begin{bmatrix}cos\theta&-sin\theta&u\\ sin\theta&cos\theta&v\\ 0&0&1\end{bmatrix}$. The mathematics behind the transformation using the matrix R us identical to step 3 in the automatic registration.

3.5 Optional ROI selection

LISA provides an optional region of interest selection after the image registration as shown in Figure 10. This step is done by connecting a polygon around the ROI to distinguish ROI from the background by setting intensity values outside of the ROI to zero. Figure 10 shows an example of optional ROI selection after registration. Comparing figure 10(a) and figure 10(b), regions at y = 7 to 15 at x = 30 of the image in figure 10(b) is trimmed off after the selection.

This step is optional to be used only when it is necessary because a user needs to be cautious when selecting the ROI at this step. If the ROI is accurately selected, it will improve the accuracy of our statistical analysis. However, if the ROI selection includes artifacts or strong noises by accident, it will corrupt the image and possibly lead to misleading significance regions.

4 Statistical Analysis

Our primary objective for LISA is to identify the areas where there have been significant changes. The statistical methods to achieve this goal are given in the following four steps, by revising the Algorithm 4.1 in our early work [13].

Curated difference 2. 2.

Non-parametric smoothing map 3. 3.

Generalized T-Map 4. 4.

FDR-controlled P-Map

4.1 Curated Difference

The initial step is to calculate curated difference of intensity values between the pairs of images in the two sequences being compared. For an image with a resolution of $n$ by $m$, the image contains $N$ number of pixels where $N$ = $n\times m$. Let $Y_{1}(i,j,t)$ and $Y_{2}(i,j,t)$ be the pixel intensity values at coordinate $(i,j)$ of the corresponding images at time $t$ in sequence 1 and sequence 2 respectively. Then the curated difference $d(i,j,t)$ is a pixel by pixel subtraction defined as $d(i,j,t)=Y(i,j,t)-Y(i,j,t)$ for $i=1,...,n$, $j=1,...,m$ and $\forall{t}$ are calculated for every pair of image comparison from the two image sequences and are used to generate the time series D-maps.

4.2 Non-parametric smoothing map

Our non-parametric smoothing map first pads the border of curated difference and then smooth the curated difference, $d_{j}$, in the region of interest. Padding the edge will prevent outstanding false-positive differences due to image alignment and noise. Our non-parametric smoothing method uses a local polynomial regression. The method is reproduced from 4.1.T-type Statistics [13] as follows.

For a single set of curated difference image data $\{(x_{1j},x_{2j},d_{j});j=1,...,N\}$, consider the model

[TABLE]

where $m(\cdot)$ is an unknown smooth function and $\epsilon_{i}$ is an error term representing the random errors of the curated differences derived from the image observations. Then the smooth function $m$ at $s=(x_{1},x_{2})$ can be approximated by the value of a fitted local polynomial function using the intensity values in the neighborhood of $s$. Specifically, for any $s=(x_{1},x_{2})$ of interest, the least squares problem seeks $\hat{\boldsymbol{\beta}}=(\hat{\beta}_{0},...,\hat{\beta}_{5})$ such that

[TABLE]

and then defines the estimate of the regression surface as $\hat{m}(x_{1},x_{2})=\hat{\beta}_{0}$ [18, 19]. The smooth weight function $w(z;h)$ peaks at $z=0$ and the smoothness of $\hat{m}$ is controlled by the bandwidth $h$. A normal density function with mean 0 and standard deviation $h$ is a typical choice. In our procedure, the bandwidths $h_{1}$ and $h_{2}$ are selected based on the approach of 10 fold cross validation with 20 iterations [20].

4.3 Generalized T-Map

The method used for generalized T-Map is reproduced from 4.1.T-type Statistics [13] as follows.

The hypothesis we are testing is $H_{0}:m(s)=0$ vs. $H_{1}:m(s)>0$ at $s=(x_{1},x_{2})$. Where $m$ can be written as a linear combination of the actual curated differences $d_{j}$ of pixel values,

[TABLE]

where $p(s)^{T}=(p_{1}(s),...,p_{N}(s))$ are the rows of the hat matrix specified by the local quadratic approximation. The estimated standard deviation of the local estimate $\hat{m}$ is $\hat{S}(s)=\hat{\sigma}\lVert(s)\rVert$. From these estimates, a T-type statistics is derived as,

[TABLE]

$T(s)$ follows t-distribution and the null hypothesis is rejected at $s$ when $T(s)>t_{1-\alpha}(\sigma_{1}^{2}/\sigma_{2})$ under a significance level of $\alpha$ with the degrees of freedom $n-6$, and $\sigma_{1}$ and $\sigma_{2}$ obtained by two-moment chi-square approximations [18, 21, 22]. T-type statistics is different from the conventional T-statistics because it uses the weighted average of the values near $s$ whereas conventional T-statistics uses simple average. Therefore, $T(s)$ and $T(s^{\prime})$ can be correlated for our T-type statistics whereas the test statistics are independent to each other in the conventional T-tests.

4.4 FDR-controlled P-Map

The P-Map is generated from multiple comparison controlled p-values derived from the distribution of the T-type statistics. The most important problem in generating P-Map is overcoming the multiplicity problem for simultaneous comparisons. Family-wise error rate (FWER) such as Bonferroni correction and false discovery rate(FDR) are the most common methods used to control for multiple comparisons. Section 4.2 of Zhang(2005)[23] discussed the relationship between FWER and FDR and we decided to use FDR for LISA. Specifically, we used the Benjamini and Hochberg’s FDR-controlled method [24]. The details of the method are introduced in section 4.2 Multiple testing problem of Wang et al.(2006)[13]. LISA provides ‘2D’ and ‘3D’ options in displaying the P-map.

5 Application

In this section, an application of LISA on Neuromuscular Electrical Stimulation(NMES) experiment [14] is introduced. NMES is a treatment proposed to improve pressure distribution at the seating or lying support area for patients with mobile difficulties. Patients are at a risk of developing pressure ulcers when seated or lying down for a long time due to the blocking of the blood flow. NMES is intended to improve the blood flow for healthier regional tissue. In this application, the effect of NMES is evaluated by comparing the sequences of the pressure mat images before and after the treatment. A step by step process of using LISA will be introduced to conduct the comparison.

Generic LISA function :

⬇

lisa(’filename1’,’filename2’,’scantype’,’nframe’,’nrow’,’ncol’,

   ’rowid’,’colid’,’preprocess’,’segtype’,’registtype’,’roi’,

   ’displaym’,’dimenpmap’,’parcomp’)

5.1 Read in Data

The argument ‘scantype’ takes one of the three possible values. When scantype = ‘row’, LISA finds the first row that matches ‘rowid’, then ‘nrow’ rows after the ‘rowid’ are stored as a frame. The search continues until the frame after the last ‘rowid’ is read in. When scantype = ‘col’, LISA scans the csv file, and read in frames using the frame numbers found under the ‘colid’ column. When scantype = ‘no’, LISA uses blank rows as row ids, and read in ‘nrow’ rows after each blank row as individual frames. After all data are read in, LISA will have two arrays with nframe number of cells having ‘nrow’ by ‘ncol’ matrix of the intensity values in each cell.

5.2 Image Processing

A user can skip the image processing step if it’s unnecessary by setting ‘preprocess’ to ‘T’. Otherwise, a user can set ‘preprocess’ to ‘F’ to process the images as follows.

5.3 Select Region of Interest

LISA shows the first images of each sequence to select ROI as the first step of image processing. An example is displayed in Figure 11 and specific instructions on how to choose ROI are provided on the screen. A user needs to drag a square to roughly select the ROI. It is advised to choose a region large enough to enclose all ROI. Between the two regions selected from each sequence, the size of the larger region will be used to crop the images.

5.4 Image Segmentation

In the image segmentation step, noises to ROI are identified using the pixel intensity values. Automatic and manual image segmentation options are available in LISA as described in section 3.3. The following is an example.

(Automatic segmentation)

Select ROI : A user can use a multi-point polygon to select detailed ROI as shown in Figure 12.

Select groups for mixture of Gaussian distribution : A user can change the number of bins and define visually distinct groups from the distribution of intensity as shown in Figure 13.

(Manual segmentation)

Pick cutoff value : A user should select one cutoff value from the histogram of the intensity values to distinguish the noise and ROI as shown in Figure 14.

5.5 Image Registration

Image registration can be done automatically or manually as described in section 3.4. Automatic registration does not require any additional input from a user. Manual registration requires a user to select two points : a reference point and a second point to form a reference line using the first image of each sequence as shown in Figure 15.

5.6 Image Alignment Check

Before conducting statistical comparison, LISA overlaps the first images of the two sequences which have been processed to confirm whether the image segmentation and image registration is satisfied by a user. If a user is satisfied and selects ‘Yes’, then LISA continues to statistical analysis but if a user selects ‘No’ LISA provides suggestions to the user on image processing in a pop-up window and terminates the process. An example of the overlapped image is shown in Figure 16.

5.7 Output Movies

Depending on the option a user selects, LISA can provide 8 different movies :

O-movie 1 (Figure 17-a) : Movie of the actual images from the first sequence 2. 2.

O-movie 2 (Figure 17-b) : Movie of the actual images from the second sequence 3. 3.

R-movie 1 (Figure 17-c) : Movie of the registered images from the first sequence 4. 4.

R-movie 2 (Figure 17-d) : Movie of the registered images from the second sequence 5. 5.

D-movie (Figure 17-e) : Movie of the curated difference 6. 6.

S-movie (Figure 17-f) : Movie of the non-parametrically smoothed difference 7. 7.

T-movie (Figure 17-g) : Movie of the T-type statistics 8. 8.

P-movie (Figure 17-h) : Movie of FDR-controlled p-values

When ‘basic’ is selected only P-map will be generated but when ‘all’ is selected, all 8 movies will be generated for display. In D-movie, S-movie and T-movie, a blue peak signifies positive differences which indicate higher pressure in the first sequence compared to the second sequence, whereas, a red peak signifies the opposite. In P-movie, red areas are statistically significantly activated/different areas. The snapshots of the movies resulting from our data application using the seated pressure mats are introduced in Figure 17.

6 Discussion

LISA is an image analysis package that processes and analyzes two sequences of images to find significant activation regions between them. In LISA, the essential components of image processing functionality include ROI selection, image segmentation and image registration and that of statistical analysis functionality include non-parametric mapping, smoothing and multiple comparison adjusted testing.

The current version of LISA was developed for and focused on sequence of images from seated pressure mats but the package can be extended to analyze other types of time series images such as satellite images, hyper-spectral chemical images, longitudinal magnetic resonance images. The interactive function during the image processing step makes it easier for users to use the package and provides users with visual evidence to derive more accurate image alignment and statistical comparison results. Statistical analysis is designed to reduce possible errors resulting from edge effects and to account for multiple comparisons using false discovery rate which makes the final result conservative and accurate. LISA also has a parallel computing function to help with computing issues for big data. High-resolution or long image sequences can be analyzed with the parallel computing function in a multi-core processing environment.

Depending on a user’s needs and preference, other image processing tools can enhance the overall analysis experience, such as the Image Processing Toolbox(IPT) within MATLAB. Among the many options in IPT, a function such as multi-point registration could help users to derive more accurate image alignments depending on the characteristics of the image. The underlying method is identical to manual registration in LISA, but users can select multi-reference points for image translation. For details, please refer to the MathWorks website [25]. Any other tools can be used for image processing as long as the images are well aligned and the results are saved in ‘csv’ format to be consumed by LISA.

Development of a standalone MATLAB package LISA as an ensemble of elements from the LASR algorithm [13] with additional features for image processing and graphical user interface is the start of our project in building an image comparison analysis package for large dynamic sequence of high-resolution images. The current version can handle relatively simple and static sequence of images. However, we aim to develop a next generation of LISA to be able to handle more complex and dynamic images, such as the entire body of a human and comparing images at many different time points. If we take a pressure map of a entire body of a human as an example, our future work will focus on segmented image processing and relevant statistical comparisons. This will allow us to compare and find significant activation regions independently separating by upper body, lower body, arms, legs and for any specific areas.

Another area of our future development includes reducing the computing time. Even with the parallel computing option, the computing time for LISA can be relatively long depending on the size of the image and sequence and the number of cores in the computer because the image comparison is done pixel-by-pixel. Pixel-by-pixel comparison is the most comprehensive and exhaustive way to compare images but can also provide most accurate results relative to using measures. Therefore, we will be aiming to develop methods such as to distribute the statistical calculation to reduce the computing time working under the essential framework of pixel-by-pixel comparison.