io.github.rocsg.fijirelax.mrialgo.MRUtils Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of fijirelax Show documentation
Show all versions of fijirelax Show documentation
FijiRelax : 3D+t MRI analysis and exploration using multi-echo spin-echo sequences
/*
*
*/
package io.github.rocsg.fijirelax.mrialgo;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Random;
import java.util.concurrent.atomic.AtomicInteger;
import org.apache.commons.math3.distribution.ChiSquaredDistribution;
import io.github.rocsg.fijiyama.registration.TransformUtils;
import io.github.rocsg.fijiyama.common.VitimageUtils;
import ij.IJ;
import ij.ImagePlus;
import ij.plugin.Duplicator;
import ij.process.ImageProcessor;
import io.github.rocsg.fijirelax.curvefit.LMDualCurveFitterNoBias;
import io.github.rocsg.fijirelax.curvefit.SimplexDualCurveFitterNoBias;
/**
* This class is a utility class to hold many helpers used for simulating MRI echoes, computing rice noise estimation, reading information into TIFF metadata.
* 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 maxMag)maxMag=echoVals[i][pix];
}
}
}
}
for(int i=0;imaxTr && Tr[i]<10000)maxTr=Tr[i];
if(Te[i]>maxTe)maxTe=Te[i];
}
tab[0]=maps[0];//PD
tab[nMapsOutput-1]=VitimageUtils.nullImage(tab[0]);//Mask
if(hasT1T2) {tab[1]=maps[1];tab[2]=maps[2];}
else if(hasT1)tab[1]=maps[1];
else if(hasT2)tab[1]=maps[2];
ImagePlus hyper=VitimageUtils.hyperStackingChannels(tab);
for(int c=nMapsOutput;c list=new ArrayList();
for(double dou : tab)list.add(dou);
Collections.sort(list);
for(int i=0;i1) {
IJ.showMessage("Critical fail in MRUtils : get "+nbCat+" categories instead of 1 \nwhen calling getDataTypeOfThisMagneticResonanceSlice(ImagePlus img,"+c+","+z+","+f+"\nLabel was : "+label);
}
if(nbCat<1)return MRDataType.OTHER;
return dataT;
}
/**
* Read values of the capillary in slices label
*
* @param img the image to test
* @param c the channel
* @param z the zslice
* @param t the timeframe
* @return the double[]
*/
public static double[] readCapValuesInSliceLabel(ImagePlus img,int c, int z, int t) {
double []ret=new double[2];
String label=img.getStack().getSliceLabel(VitimageUtils.getCorrespondingSliceInHyperImage(img, c, z, t));
String[]strTab=label.split("_");
for(String str:strTab) {
if(str.contains("CAP")) {
ret[0]=Double.parseDouble((str.split("=")[1]).split("\\+-")[0]);
ret[1]=Double.parseDouble((str.split("=")[1]).split("\\+-")[1]);
}
}
return ret;
}
/**
* Compute the mean of the double tab
*
* @param tab the tab
* @return the double
*/
protected static double mean(double[] tab) {
double m=0;
for(int i=0;i1) {
if(prm0Sigma_1Tr_2Te==PARAM_SIGMA && prm[0].equals("SIGMARICE"))return Double.parseDouble(prm[1]);
if(prm0Sigma_1Tr_2Te==PARAM_TR && prm[0].equals("TR"))return Double.parseDouble(prm[1]);
if(prm0Sigma_1Tr_2Te==PARAM_TE && prm[0].equals("TE"))return Double.parseDouble(prm[1]);
}
}
return 0;
}
/**
* Modify Rice noise sigma, and write it in the slice label
*
* @param hyperImg the hyper img
* @param c the channel
* @param z the zslice
* @param t the timeframe
* @param newSigma the new sigma to write
*/
public static void modifySigma(ImagePlus hyperImg,int c,int z,int t,double newSigma) {
String label=hyperImg.getStack().getSliceLabel(VitimageUtils.getCorrespondingSliceInHyperImage(hyperImg, c, z, t));
String[]strTab=label.split("_");
String newLab="";
for(int i=0;itab[i]*1.3 || tabTemp[j]xMax-1) || (yCor>yMax-1)) {IJ.log("Bad coordinates. Data set to 0"); return null;}
for(int z=zm;z<=zM;z++) {
for(int x=xm;x<=xM;x++) {
for(int y=ym;y<=yM;y++) {
for(int ec=0;ecxMax-1) || (yCor>yMax-1)) {IJ.log("Bad coordinates. Data set to 0"); return null;}
for(int z=zm;z<=zM;z++) {
for(int x=xm;x<=xM;x++) {
for(int y=ym;y<=yM;y++) {
for(int ec=0;ecxMax-1) || (yCor>yMax-1)) {IJ.log("Bad coordinates. Data set to 0"); return data;}
double[][][]weights=new double[xM-xm+1][yM-ym+1][zM-zm+1];
double sigmaX=crossWidth;
double sigmaY=crossWidth;
double sigmaZ=crossThick;
double sum=0;
for(int x=xm;x<=xM;x++) {
for(int y=ym;y<=yM;y++) {
for(int z=zm;z<=zM;z++) {
if(!gaussianWeighting)weights[x-xm][y-ym][z-zm]=1.0/nHits;
else {
if(sigmaX<=0 && sigmaZ<=0) {//One point
weights[x-xm][y-ym][z-zm]=1;
}
else if(sigmaX<=0 && sigmaZ>0) {//Mono dimensional along z
weights[x-xm][y-ym][z-zm]=1/(Math.pow(2*Math.PI,0.5)*sigmaZ) * Math.exp(- (z-zCor)*(z-zCor)/(sigmaZ*sigmaZ));
}
else if(sigmaX>0 && sigmaZ<=0) {//Two dimensional along x and y
weights[x-xm][y-ym][z-zm]=1/(Math.pow(2*Math.PI,1)*sigmaX*sigmaY) * Math.exp(- (x-xCor)*(x-xCor)/(sigmaX*sigmaX) - (y-yCor)*(y-yCor)/(sigmaY*sigmaY));
}
else {//Three dimensional gaussian
weights[x-xm][y-ym][z-zm]=1/(Math.pow(2*Math.PI,1.5)*sigmaX*sigmaY*sigmaZ) * Math.exp(- (x-xCor)*(x-xCor)/(sigmaX*sigmaX) - (y-yCor)*(y-yCor)/(sigmaY*sigmaY) - (z-zCor)*(z-zCor)/(sigmaZ*sigmaZ));
}
}
sum+=weights[x-xm][y-ym][z-zm];
}
}
}
if(gaussianWeighting) {
for(int x=xm;x<=xM;x++) {
for(int y=ym;y<=yM;y++) {
for(int z=zm;z<=zM;z++) {
weights[x-xm][y-ym][z-zm]/=sum;
}
}
}
}
for(int ec=0;ec0)for(int p=0;p9999 ? param[3] : 0)) / param[2]) ),sigma);break;
case T1T2_DEFAULT_T2_MULTI_RICE:ret=RiceEstimator.besFunkCost( (param[0]*Math.exp(-( (techo+(trecov>9999 ? param[5] : 0)) / param[3]) )+param[1]*Math.exp(-((techo+(trecov>9999 ? param[5] : 0)) / param[4])) * (1-Math.exp(-((trecov) / param[2])) ) ) , sigma);break;
default:IJ.showMessage("In MRUtils.getFitFunctionValue, unexpected fit type : "+fitType+" . \nBrutal stop now.");ret=1/0;
}
return ret;
}
/**
* Compute normalised khi2 and p values according to the sigma parameter of the Rice noise
* the number of points averaged for the estimation is used to estimate sigma of the sum
* Another interesting point is that sigma is reestimated locally for each magnitude value
* as the rice sigma is not equal to the observed sigma for different values of the magnitude.
*
* @param tabData the tab data
* @param tabTrTimes the tab of recovery times
* @param tabTeTimes the tab of echo times
* @param sigma the sigma of Rice noise
* @param estimatedParams the estimated parameters of the fitting function
* @param fitType the fit type (ex: T2_MONO_BIAS)
* @param debug the debug flag in order to display some debugging infformations
* @param riceEstimator the RiceEstimator object in use to handle the analysis
* @param nbPoints the nb of points averaged for getting this data
* @return the double[]{khi2, p}
*/
public static double[] fittingAccuracies(double[] tabData,double[] tabTrTimes,double[] tabTeTimes,double sigma,double[] estimatedParams,int fitType,boolean debug,RiceEstimator riceEstimator,int nbPoints){
int N_PARAMS=getNparams(fitType);
double []realP=fittenRelaxationCurve(tabTrTimes,tabTeTimes,estimatedParams,sigma,fitType);
double khi2=0; double diff=0; double sigmaEst=0;
for(int indT=0;indT(X-radiusAroundCapillary-5) || coords[z][0]>(Y-radiusAroundCapillary-5) ) {
IJ.showMessage("Warning in normalizeBeforeComputation : unsteady mean computation at z="+z+" because capillary is detected near the border, at coordinates "+coords[z][0]+" , "+coords[z][1]+" with an indicated security radius of "+radiusAroundCapillary);
}
//Mesurer la moyenne autour du capillaire
meanOverCap[z]=VitimageUtils.meanValueofImageAround(sliceTemp, coords[z][0], coords[z][1], 0,radiusAroundCapillary);
//Mesurer la moyenne de l'image
int radiusImg=Math.min(X/2,Y/2)-2;
meanOverAll[z]=VitimageUtils.meanValueofImageAround(sliceTemp, X/2, Y/2, 0,radiusImg);
//Oter l'un à l'autre
meanOverRest[z]=meanOverAll[z]*(2*radiusImg+1)*(2*radiusImg+1) - meanOverCap[z]*(2*radiusAroundCapillary+1)*(2*radiusAroundCapillary+1);
meanOverRest[z]=meanOverRest[z]/( (2*radiusImg+1)*(2*radiusImg+1) - (2*radiusAroundCapillary+1)*(2*radiusAroundCapillary+1) );
}
//Faire la moyenne des valeurs mesurées along Z
globalMeanOverRest=VitimageUtils.statistics1D(meanOverRest)[0];
globalMeanOverCap=VitimageUtils.statistics1D(meanOverCap)[0];
int radiusSq=radiusAroundCapillary*radiusAroundCapillary;
double distSq;
//Pour chaque Z
for(int im=0;im