HalfValueLayerFunction.java

package edu.stanford.rsl.conrad.physics;

import edu.stanford.rsl.conrad.utils.CONRAD;
import edu.stanford.rsl.conrad.utils.LinearInterpolatingDoubleArray;
import edu.stanford.rsl.jpop.SimpleFunctionOptimizer;
import edu.stanford.rsl.jpop.SimpleOptimizableFunction;


public class HalfValueLayerFunction implements SimpleOptimizableFunction {

	private double [] xRayEnergy;
	private double [] xRaySpectrum;	
	private LinearInterpolatingDoubleArray attenuationCoefficients1;
	private LinearInterpolatingDoubleArray attenuationCoefficients2;
	private LinearInterpolatingDoubleArray attenuationCoefficients3;
	public double optimalX;
	
	public HalfValueLayerFunction (double [] energies, double [] spectrum, 
			LinearInterpolatingDoubleArray lida1, 
			LinearInterpolatingDoubleArray lida2, 
			LinearInterpolatingDoubleArray lida3){
		this.xRayEnergy = energies; // energy KeV
		this.xRaySpectrum = spectrum; // photon count #
		this.attenuationCoefficients1 = lida1;
		this.attenuationCoefficients2 = lida2;
		this.attenuationCoefficients3 = lida3;
	}

	@Override
	public double evaluate(double x) {
		
		double sumN0 = 0;
		double sumN1 = 0;
		double coeff1; // Al
		double coeff2; // CsI converter
		double coeff3; // Air absorption
		double detectorEfficiency = 1;
		double xRayQuanta = 1;
		boolean isHVLMeasuredAtDetector = false;
		boolean isAirAbsorptionCounted = true;
		
		try {
			for (int i = 0; i< xRaySpectrum.length; i++){
				coeff1 = attenuationCoefficients1.getValue(xRayEnergy[i]/1000);
				coeff2 = attenuationCoefficients2.getValue(xRayEnergy[i]/1000);
				coeff3 = attenuationCoefficients3.getValue(xRayEnergy[i]/1000);
				
				if(isHVLMeasuredAtDetector) {
					// Detector flat panel efficiency, thickness = 200 [g/cm2]
					detectorEfficiency = 1 - Math.exp(-coeff2 * 200);						
				}
				
				if(isAirAbsorptionCounted) {
					xRayQuanta = XRaySpectrum.xrqpRe(coeff3);
				}				
				
				// Exposure (R) of the spectrum specified by the vector phi at energies specified by E (keV).
				sumN1 += xRaySpectrum[i] * Math.exp(-coeff1 * x) * detectorEfficiency * xRayEnergy[i] / xRayQuanta;
				sumN0 += xRaySpectrum[i] * detectorEfficiency * xRayEnergy[i] / xRayQuanta;		
				
				// System.out.println(xRayEnergy[i]+"\t"+coeff1+"\t"+coeff2+"\t"+coeff3);
			}	
		}catch (Exception e){
			e.printStackTrace();
		}
		
		// System.out.println(x +"\t"+ Math.abs(sumN0/2 - sumN1));	
		
		return Math.abs(sumN0/2 - sumN1);
	}
	
	public void runOptimization() {
		SimpleFunctionOptimizer optimizer = new SimpleFunctionOptimizer();
		optimizer.setLeftEndPoint(0);
		optimizer.setRightEndPoint(10);
		//optimizer.setInitialGuess(0.44);
		optimizer.setAbsoluteTolerance(CONRAD.DOUBLE_EPSILON);
		optimizer.setRelativeTolerance(CONRAD.DOUBLE_EPSILON);
		//optimizer.
		optimalX = optimizer.minimize(this);		
	}
	
	public double getOptimalX(){
		return optimalX;
	}

}
/*
 * Copyright (C) 2010-2014 Andreas Maier
 * CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/