//****************************************************************************** // // File: SpinSignalDifference.java // Package: edu.rit.mri // Unit: Class edu.rit.mri.SpinSignalDifference // // 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.numeric.Series; import edu.rit.numeric.VectorFunction; import edu.rit.numeric.XYSeries; /** * Class SpinSignalDifference provides an object that computes the difference * between a measured spin signal and a model spin signal. Class * SpinSignalDifference implements interface {@linkplain * edu.rit.numeric.VectorFunction VectorFunction}. *

* An instance of class SpinSignalDifference is constructed with two data * series, each an instance of class {@linkplain edu.rit.numeric.Series Series}. * The first data series contains the measurement times * ti. The second data series contains the measured * spin signals S(ti). The vector function's * result vector length, M, is the same as the length of the data series. *

* The vector function's argument vector x gives parameters for the model * spin signal. The argument vector consists of one or more pairs of consecutive * values. Each pair of consecutive values corresponds to one tissue. The first * value of the pair is ρj, the spin density for * tissue j. The second value of the pair is Rj, * the spin relaxation rate for tissue j. The number of tissues is * specified as a constructor parameter. The vector function's argument vector * length, N, is twice the number of tissues. *

* The vector function is calculated as follows. For each i, 0 <= * i <= M−1: *

* fi(x) = Model S(ti;x) − Measured S(ti) *
*

* The model spin signal is defined as follows: *

* Model S(ti;x) = Σj ρj [1 − 2 exp(−Rj ti)] *
*

* To find the model spin signal parameters x that best fit the measured * spin signal, use an instance of class SpinSignalDifference with the nonlinear * least squares algorithm in class {@linkplain * edu.rit.numeric.NonLinearLeastSquares NonLinearLeastSquares}. * * @author Alan Kaminsky * @version 11-Jun-2008 */ public class SpinSignalDifference implements VectorFunction { // Hidden data members. private int M; private int L; private Series t_series; private Series S_series; // Exported constructors. /** * Construct a new spin signal difference function. * * @param data_series X-Y series of measured time values (X) and measured * spin signal values (Y). * @param L Number of tissues in the model spin signal. * * @exception NullPointerException * (unchecked exception) Thrown if data_series is null. * @exception IllegalArgumentException * (unchecked exception) Thrown if the series' length is 0. Thrown if * L <= 0. */ public SpinSignalDifference (XYSeries data_series, int L) { this (data_series.xSeries(), data_series.ySeries(), L); } /** * Construct a new spin signal difference function. * * @param t_series Series of measured time values. * @param S_series Series of measured spin signal values. * @param L Number of tissues in the model spin signal. * * @exception NullPointerException * (unchecked exception) Thrown if t_series is null or * S_series is null. * @exception IllegalArgumentException * (unchecked exception) Thrown if either series' length is 0. Thrown if * the two series have different lengths. Thrown if L <= 0. */ public SpinSignalDifference (Series t_series, Series S_series, int L) { if (t_series == null) { throw new NullPointerException ("SpinSignalDifference(): t_series is null"); } if (S_series == null) { throw new NullPointerException ("SpinSignalDifference(): S_series is null"); } if (t_series.length() == 0) { throw new IllegalArgumentException ("SpinSignalDifference(): t_series is zero length, illegal"); } if (S_series.length() == 0) { throw new IllegalArgumentException ("SpinSignalDifference(): S_series is zero length, illegal"); } if (t_series.length() != S_series.length()) { throw new IllegalArgumentException ("SpinSignalDifference(): t_series and S_series have different lengths, illegal"); } if (L <= 0) { throw new IllegalArgumentException ("SpinSignalDifference(): L (= "+L+") <= 0, illegal"); } this.M = t_series.length(); this.L = L; this.t_series = t_series; this.S_series = S_series; } // Exported operations. /** * Returns the length of the result vector, M. * * @return M. */ public int resultLength() { return M; } /** * Returns the length of the argument vector, N. * * @return N. */ public int argumentLength() { return L<<1; } /** * Evaluate this function with the given argument vector. The result is * stored in the vector y. Specifically, for i = 0 to * M−1, yi = * fi(x). * * @param x Argument vector (input). Must be an N-element array. * @param y Result vector (output). Must be an M-element array. * * @exception DomainException * (unchecked exception) Thrown if any argument in x is outside * the allowed set of values for this function. * @exception RangeException * (unchecked exception) Thrown if any element of the result vector is * outside the range of type double. */ public void f (double[] x, double[] y) { // Subtract out the measured spin signal. for (int i = 0; i < M; ++ i) { y[i] = -S_series.x(i); } // Add in the model spin signal. for (int i = 0; i < M; ++ i) { double t_i = t_series.x(i); for (int j = 0; j < L; ++ j) { double dens_j = x[(j<<1)]; double rate_j = x[(j<<1)+1]; y[i] += SpinSignal.S (dens_j, rate_j, t_i); } } } /** * Calculate this function's Jacobian matrix with the given argument vector. * The result is stored in the matrix dydx. Specifically, for * i = 0 to M−1 and j = 0 to N−1, * dydxi,j = * ∂fi(x) ⁄ ∂xj . * * @param x Argument vector (input). Must be an N-element array. * @param dydx Jacobian matrix (output). Must be an * M×N-element matrix. * * @exception DomainException * (unchecked exception) Thrown if any argument in x is outside * the allowed set of values for this function. * @exception RangeException * (unchecked exception) Thrown if any element of the Jacobian matrix is * outside the range of type double. */ public void df (double[] x, double[][] dydx) { for (int i = 0; i < M; ++ i) { double t_i = t_series.x(i); double[] dydx_i = dydx[i]; for (int j = 0; j < L; ++ j) { double dens_j = x[(j<<1)]; double rate_j = x[(j<<1)+1]; double twoexp = 2.0*Math.exp(-rate_j*t_i); dydx_i[(j<<1)] = 1.0 - twoexp; dydx_i[(j<<1)+1] = dens_j*t_i*twoexp; } } } }