//******************************************************************************
//
// File: PixelAnalysis.java
// Package: edu.rit.mri
// Unit: Class edu.rit.mri.PixelAnalysis
//
// This Java source file is copyright (C) 2008 by Alan Kaminsky. All rights
// reserved. For further information, contact the author, Alan Kaminsky, at
// ark@cs.rit.edu.
//
// This Java source file is part of the Parallel Java Library ("PJ"). PJ is free
// software; you can redistribute it and/or modify it under the terms of the GNU
// General Public License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// PJ is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
// A PARTICULAR PURPOSE. See the GNU General Public License for more details.
//
// Linking this library statically or dynamically with other modules is making a
// combined work based on this library. Thus, the terms and conditions of the
// GNU General Public License cover the whole combination.
//
// As a special exception, the copyright holders of this library give you
// permission to link this library with independent modules to produce an
// executable, regardless of the license terms of these independent modules, and
// to copy and distribute the resulting executable under terms of your choice,
// provided that you also meet, for each linked independent module, the terms
// and conditions of the license of that module. An independent module is a
// module which is not derived from or based on this library. If you modify this
// library, you may extend this exception to your version of the library, but
// you are not obligated to do so. If you do not wish to do so, delete this
// exception statement from your version.
//
// A copy of the GNU General Public License is provided in the file gpl.txt. You
// may also obtain a copy of the GNU General Public License on the World Wide
// Web at http://www.gnu.org/licenses/gpl.html.
//
//******************************************************************************
package edu.rit.mri;
import edu.rit.mri.SpinSignal;
import edu.rit.mri.SpinSignalDifference;
import edu.rit.numeric.NonLinearLeastSquares;
import edu.rit.numeric.NonNegativeLeastSquares;
import edu.rit.numeric.Series;
import edu.rit.numeric.TooManyIterationsException;
import java.util.ArrayList;
import java.util.List;
/**
* Class PixelAnalysis provides a routine for doing a spin relaxometry analysis
* on one pixel of a magnetic resonance image.
*
* The input to the analysis is a measured spin signal expressed as two
* {@linkplain edu.rit.numeric.Series Series} objects, a time series
* ti and a spin signal series
* S(ti). Another input is a series of fixed
* spin-lattice relaxation rates R1j. These rates are
* chosen to cover the range of likely rates for the magnetic resonance image
* being analyzed.
*
* The routine first does a nonnegative, linear least squares fit of the
* input data to a model consisting of a group of tissues with the input fixed
* spin-lattice relaxation rates R1j. Peaks in the
* linear least squares fit determine the number of tissues and the approximate
* spin density and spin-lattice relaxation rate for each tissue.
*
* The routine then does a nonlinear least squares fit of the input data
* to a model consisting of the number of tissues determined in the previous
* step. The nonlinear least squares fit "polishes up" the spin densities and
* spin-lattice relaxation rates determined in the previous step, which are only
* approximate.
*
* The routine checks the nonlinear least squares fit for plausibility. To be
* plausible:
*
* -
* All spin densities and spin-lattice relaxation rates must be positive. (If
* this is not the case, it's likely the nonlinear least squares fit is trying
* to fit the data to too many parameters, resulting in nonsensical parameter
* values.)
*
-
* All the spin-lattice relaxation rates must be sufficiently far apart.
* Specifically, the relative difference between any two rates must be greater
* than 0.001. (If the rates are closer together than that, it's likely they
* represent the same tissue.)
*
-
* The sum of the spin densities must agree with the asymptotic spin signal for
* large values of t. Specifically, the sum must be within 20% of the
* average of the last seven spin signal values. (If this is not the case, it's
* likely that two spurious "tissues," one with a very large relaxation rate,
* one with a very small relaxation rate, are canceling each other out, and
* there really should be only one tissue.)
*
* If the nonlinear least squares fit is not plausible, the routine decides it
* is trying to fit too many tissues. The routine eliminates the tissue with the
* smallest spin density and repeats the nonlinear least squares fit. This
* continues until the fit is plausible or until all the tissues have been
* eliminated, in which case the routine reports that it could not find a
* solution.
*
* The output of the analysis is a list of the tissues' computed spin densities
* and a list of the tissues' computed spin-lattice relaxation rates.
*
* @author Alan Kaminsky
* @version 16-Jun-2008
*/
public class PixelAnalysis
{
// Prevent construction.
private PixelAnalysis()
{
}
// Exported operations.
/**
* Do a spin relaxometry analysis.
*
* @param t_series
* Series of measured time values, of length M (input).
* @param S_series
* Series of measured spin signal values, of length M (input).
* @param R1_series
* Series of fixed spin-lattice relaxation rates for the linear part of
* the analysis, of length N (input).
* @param A
* Design matrix for the linear part of the analysis (input). This must
* be an M×N-element matrix such that
* Ai,j = 1 − 2
* exp(−R1j ti).
* (The design matrix is supplied as an argument because the same design
* matrix is typically used for every pixel in an image, and calculating
* the design matrix just once outside this routine saves time.)
* @param rho_list
* List in which to store the computed spin densities (output). The size
* of the list is the number of tissues. If the routine could not find a
* solution, the size of the list is 0.
* @param R1_list
* List in which to store the computed spin-lattice relaxation rates
* (output). The size of the list is the number of tissues. If the
* routine could not find a solution, the size of the list is 0.
*/
public static void analyze
(Series t_series,
Series S_series,
Series R1_series,
double[][] A,
List rho_list,
List R1_list)
{
int M = t_series.length();
int N = R1_series.length();
// Do a spin relaxometry analysis using nonnegative linear least
// squares.
// Create nonnegative linear least squares solver.
NonNegativeLeastSquares linsolver = new NonNegativeLeastSquares (M, N);
// Find the solution.
for (int i = 0; i < M; ++ i)
{
System.arraycopy (A[i], 0, linsolver.a[i], 0, N);
linsolver.b[i] = S_series.x(i);
}
linsolver.solve();
double[] rho_series = linsolver.x;
// Find peaks in the solution. A peak occurs at index i if
// rho[i] > rho[i-1] and rho[i] > rho[i+1].
ArrayList approx_rho_list = new ArrayList();
ArrayList approx_R1_list = new ArrayList();
for (int j = 0; j < N; ++ j)
{
if (rho_series[j] > (j == 0 ? 0.0 : rho_series[j-1]) &&
rho_series[j] > (j == N-1 ? 0.0 : rho_series[j+1]))
{
approx_rho_list.add (rho_series[j]);
approx_R1_list.add (R1_series.x(j));
}
}
// Do a spin relaxometry analysis using nonlinear least squares. Peaks
// in the linear analysis give the initial vector of densities and
// rates.
// Repeat until the solution is plausible.
boolean plausible = false;
int L = approx_rho_list.size();
rho_list.clear();
R1_list.clear();
while (L > 0 && ! plausible)
{
// Create spin signal difference function. L = number of tissues.
SpinSignalDifference fcn =
new SpinSignalDifference (t_series, S_series, L);
// Create nonlinear least squares solver.
NonLinearLeastSquares nonlinsolver =
new NonLinearLeastSquares (fcn);
// Find the solution.
for (int i = 0; i < L; ++ i)
{
nonlinsolver.x[(i<<1)] = approx_rho_list.get(i);
nonlinsolver.x[(i<<1)+1] = approx_R1_list.get(i);
}
try
{
nonlinsolver.solve();
for (int i = 0; i < L; ++ i)
{
rho_list.add (nonlinsolver.x[(i<<1)]);
R1_list.add (nonlinsolver.x[(i<<1)+1]);
}
// Decide if solution is plausible.
plausible = checkPlausibility (S_series, rho_list, R1_list);
}
// Couldn't find a solution.
catch (TooManyIterationsException exc)
{
plausible = false;
}
// If solution is not plausible, eliminate tissue with smallest
// density and try again.
if (! plausible)
{
double minrho = Double.MAX_VALUE;
int mini = 0;
for (int i = 0; i < L; ++ i)
{
if (approx_rho_list.get(i) < minrho)
{
minrho = approx_rho_list.get(i);
mini = i;
}
}
approx_rho_list.remove (mini);
approx_R1_list.remove (mini);
L = approx_rho_list.size();
rho_list.clear();
R1_list.clear();
}
}
}
// Hidden operations.
/**
* Decide if the given solution is plausible.
*/
private static boolean checkPlausibility
(Series S_series,
List rho_list,
List R1_list)
{
int M = S_series.length();
int L = rho_list.size();
// If any density or rate is negative, solution is not plausible.
for (int i = 0; i < L; ++ i)
{
if (rho_list.get(i) < 0.0)
{
return false;
}
if (R1_list.get(i) < 0.0)
{
return false;
}
}
// If relative difference between any two rates is too small, solution
// is not plausible.
for (int i = 0; i < L-1; ++ i)
{
double R_i = R1_list.get(i);
for (int j = i+1; j < L; ++ j)
{
double R_j = R1_list.get(j);
double reldiff = 2.0*Math.abs(R_i-R_j)/Math.abs(R_i+R_j);
if (reldiff <= 0.001)
{
return false;
}
}
}
// If sum of densities is too far from asymptotic measurement for large
// t, solution is not plausible.
double sumrho = 0.0;
for (int i = 0; i < L; ++ i)
{
sumrho += rho_list.get(i);
}
double S_last = 0.0;
int n = 0;
for (int i = M-1; i >=0 && n < 7; -- i)
{
S_last += S_series.x(i);
++ n;
}
S_last /= n;
double reldiff = Math.abs(sumrho-S_last)/Math.abs(S_last);
if (reldiff >= 0.2)
{
return false;
}
// Solution is plausible.
return true;
}
}