• Home
• io.github.rocsg
• fijirelax
• 4.0.4
• source code
• RiceEstimator.java

×

Project download

Many resources are needed to download a project. Please understand that we have to compensate our server costs. Thank you in advance.
Project price only 1 $

You can buy this project and download/modify it how often you want.

io.github.rocsg.fijirelax.mrialgo.RiceEstimator Maven / Gradle / Ivy

package io.github.rocsg.fijirelax.mrialgo;

import java.util.Random;

import io.github.rocsg.fijiyama.common.VitimageUtils;

/** This class provides utilities to estimate MRI relaxation parameters in presence of a Rice noise
 * Rice noise affects the signal in a way that make it complicated to invert : its moments depend on the value of the unaltered signal
 * For more information, refer to Fernandez et al. 2022 FijiRelax: Fast and noise-corrected estimation of MRI relaxation maps in 3D + t (in prep.)
 * 
 * 
 * 	Copyright (C) 2022  
 *
 *   This program 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.
 *
 *   This program 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.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see 
*/
public class RiceEstimator {


	/** The Constant epsilon. */
	private final static double epsilon=0.00001;
	
	/** The Constant debug. */
	private final static boolean debug=false;
	
	/** The sigma rice min. */
	private double sigmaRiceMin;
	
	/** The sigma rice max. */
	private double sigmaRiceMax;
	
	/** The observation min. */
	private double observationMin;
	
	/** The observation max. */
	private double observationMax;
	
	/** The observation range. */
	private double []observationRange;
	
	/** The sigma range. */
	private double []sigmaRange;
	
	/** The n sig. */
	private int nSig;
	
	/** The n obs. */
	private int nObs;
	
	/** The lut. */
	private double[][][]lut;
	
	

	
	/**
	 * Test of the random Rice simulator
	 */
	public static double[] testRandomRiceMeanDrift() {
		int N=100000;
		int SNR=2;
		double signal=700;
		double sigma=signal/SNR;
		double[]tab=new double[N];
		for(int i=0;i=nSig-1)coord=nSig-1-epsilon;
		return coord;
	}
	
	/**
	 * Builds the observation range.
	 */
	public void buildObservationRange() {
		if(debug)System.out.print("Construction array observations : ");
		double multFactor=this.observationMax/this.observationMin;
		for(int obs=0;obs=this.observationMax)return nObs-1-epsilon;
		else if(observation<=this.observationMin)return epsilon;
		
		double valMin=this.observationMin;double valMax=this.observationMax;int indMin=0;int indMax=nObs-1;int indMed;double valMed;
		do {
			indMed=(int)Math.round(0.5*(indMax+indMin));
			valMed=observationRange[indMed];
			if(observation1);
		return indMin+(observation-valMin)/(valMax-valMin);
	}

	
	/**
	 * Builds the lookup table among possible values of sigma and observations
	 */
	public void buildLookupTable() {
		double interMin, interMax,interMed,valMin,valMax,valMed;
		System.out.print("Building lookup table for estimation of initial signal amplitude using max likelyhood reverse estimation from observations X rice noise ...");
		int nTot=this.nObs*this.nSig;
		int incr=0;
		this.lut=new double[nSig][nObs][3];//SigmaRice, valObservation, EquivalentNormalSigma
		for(int sig=0;sigcurObs) {}
				else {
					//Recherche iterative
					valMin=besFunkCost(interMin, curSigRice);
					valMax=besFunkCost(interMax, curSigRice);
					do {
						interMed=(interMax+interMin)/2;
						valMed=besFunkCost(interMed, curSigRice);
						if(valMed>=curObs) {interMax=interMed;valMax=valMed;}
						else {interMin=interMed;valMin=valMed;}
					}
					while(Math.abs(valMed-curObs)>epsilon);
				}
				lut[sig][obs]=new double[] {curSigRice,curObs,interMed};
				if(debug)System.out.println("sig"+sig+"  obs"+obs+"  : sigma="+curSigRice+"  observation="+ curObs+"  initialSignal="+interMed);
			}
		}
		System.out.println();
	}
			
	
	/**
	 * Estimate the besselFunction sigma for these values
	 *
	 * @param sigmaRice the sigma rice
	 * @param observation the observation
	 * @return the double
	 */
	public double estimateSigma(double sigmaRice,double observation) {
		double initialSignal=estimateInitialSignal(sigmaRice,observation);
		return besFunkSigma(initialSignal,sigmaRice);
	}
	
	/**
	 * Estimate the besselFunction sigma for these values
	 *
	 * @param sigmaRice the sigma rice
	 * @param observations the observations
	 * @return the double[]
	 */
	public double []estimateSigmas(double sigmaRice,double []observations) {
		double[]ret=new double[observations.length];
		for(int r=0;r0.3) System.out.println("Warning : Acquisition > computeRiceSigmaStatic. Given :M="+meanBg+" , S="+sigmaBg+" gives sm="+val1+" , ss="+val2+"  .sm/ss="+VitimageUtils.dou(val1/val2)+". Using the first one..");
		
		return val1;
	}

	

	/**
	 * Estimation of the rice noise stddev from values measured in the background
	 * assuming that there is no object in such background. The value is computed by two methods, using the mean and stdev of the BG, then the two values are returned
	 *
	 * @param sigmaRice the sigma rice
	 * @return A double array yielding the estimated sigma in the two methods
	 */
	public static double []computeSigmaAndMeanBgFromRiceSigmaStatic(double sigmaRice) {
		boolean debug=true;
		double meanBg=sigmaRice/(Math.sqrt(2.0/Math.PI));
		double stdBg=sigmaRice/(Math.sqrt(2.0/(4.0-Math.PI)));
		return new double[] {meanBg,stdBg};
	}

	
	/**
	 * Simulation of Rice Noise over an unaltered signal.
	 *
	 * @param originalSignal the original signal
	 * @param sigmaRice the sigma rice
	 * @param nRepets the n repets
	 * @return A double corresponding to the value of the altered signal
	*/ 
	public static double getRandomRiceRealization(double originalSignal,double sigmaRice,int nRepets){
		double acc=0;
		for(int i=0;i