/*
* Copyright (C) 2010-2014 Andreas Maier, Rotimi X Ojo
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/
package edu.stanford.rsl.conrad.physics.absorption;
import java.util.ArrayList;
import java.util.HashMap;
import edu.stanford.rsl.conrad.geometry.shapes.simple.Edge;
import edu.stanford.rsl.conrad.physics.PhysicalObject;
import edu.stanford.rsl.conrad.physics.PolychromaticXRaySpectrum;
import edu.stanford.rsl.conrad.physics.materials.Material;
import edu.stanford.rsl.conrad.physics.materials.utils.AttenuationType;
import edu.stanford.rsl.conrad.physics.materials.utils.MaterialUtils;
import edu.stanford.rsl.conrad.utils.StatisticsUtil;
/**
* <p>Polychromatic Absorption Model with dynamic spectrum support. <br> This class models the absorption of an X-Ray spectrum by well-defined materials. </p>
*
* @author Rotimi X Ojo, Andreas Maier
*/
public class PolychromaticAbsorptionModel extends AbsorptionModel {
/**
*
*/
private static final long serialVersionUID = -3430401648474324638L;
private AttenuationType att = AttenuationType.TOTAL_WITH_COHERENT_ATTENUATION;
protected PolychromaticXRaySpectrum inputSpectrum;
protected double [] energies;
protected double [] photonFlux;
protected HashMap<Material, double []> attenuationCoefficientsMap;
/**
* Compute the absorption of a given X-Ray spectrum along a segmented path.
* @param segments are material dependent segmentation of the path followed by an arbitrary X-Ray
* @return the log attenuation of line integral of the product of X-Ray Intensities over a range of energies and energy dependent attenuation.
*/
@Override
public double evaluateLineIntegral(ArrayList<PhysicalObject> segments) {
double intensity = computeIntensity(segments, energies[0], energies[energies.length-1], false, true);
// we do not require the pixel area here, as we normalize with the total spectral intensity.
double value = intensity/inputSpectrum.getTotalIntensity();
//Due to the accumulation round off error, it is possible for output intensity to be marginally greater
//than input intensity in vacuum. We deal with this problem by bounding output intensity by the input
//intensity.
if(value > 1){
value = 1;
//throw new RuntimeException("PolychromaticAbsorptionModel: numerical instability found.");
}
//End
return -Math.log(value);
}
private int convertToIndex(double keV){
return (int)Math.round((keV-energies[0])/(energies[1]-energies[0]));
}
/**
* Change the X-Ray spectrum used by the absorption model. This method is important for implementing Dynamic Spectrum and XML Serialization.
* @param spectrum is the new spectrum of the absorption model.
*/
public void setInputSpectrum(PolychromaticXRaySpectrum spectrum){
this.inputSpectrum = spectrum;
// precompute energies:
energies = inputSpectrum.getPhotonEnergies();
// precompute intensities:
photonFlux = new double[energies.length];
for (int i = 0; i < energies.length; i++){
photonFlux[i]=inputSpectrum.getPhotonFlux(energies[i]);
}
// precompute absorption spectra for all known materials in database:
attenuationCoefficientsMap = MaterialUtils.loadAttenuationCoefficients(energies, att);
}
/**
* Retrieve the current input X-Ray spectrum used by the absorption model
* @return the current X-Ray spectrum used by the absorption model
*/
public PolychromaticXRaySpectrum getInputSpectrum(){
return inputSpectrum;
}
@Override
public String toString() {
return "Polychromatic Absorption Model";
}
public double getMaximalEnergy(){
return energies[energies.length-1];
}
public double getMinimalEnergy(){
return energies[0];
}
@Override
public void configure() throws Exception {
inputSpectrum = new PolychromaticXRaySpectrum(2.5);
// we use this method to preload different material properties in order to be able to execute the code in parallel efficiently.
setInputSpectrum(inputSpectrum);
}
@Override
public boolean isConfigured() {
return true;
}
public double getPhotonFluxIntegral(double startEnergy, double stopEnergy) {
int start = convertToIndex(startEnergy);
int end = convertToIndex(stopEnergy);
double sum = 0;
if (start < 0) start =0;
if (end >= photonFlux.length) end = photonFlux.length-1;
for (int i= start; i<=end; i++){
sum+=this.photonFlux[i];
}
return sum;
}
/**
* Computes the integral over the spectrum given a start and an end energy.
* If energyIntegrating is true the accumulated energy is returned otherwise the photon count.
* @param segments the path segments
* @param startEnergy the start energy [keV]
* @param endEnergy the end energy [keV]
* @param noise if true noise is generated
* @param energyIntegrating if true the photon count is weighted with the energy
* @return the intensity.
*/
public double computeIntensity(ArrayList<PhysicalObject> segments,
double startEnergy, double endEnergy, boolean noise, boolean energyIntegrating) {
double intensity = 0;
int start = convertToIndex(startEnergy);
int end = convertToIndex(endEnergy);
if (start < 0) start =0;
if (end >= photonFlux.length) end = photonFlux.length-1;
double lens [] = new double[segments.size()];
for (int j = 0; j < segments.size(); j++){
lens[j] = ((Edge)segments.get(j).getShape()).getLength();
}
for (int e = start; e <= end; e++){
double sum = 0;
for (int j = 0; j < segments.size(); j++){
PhysicalObject o = segments.get(j);
double att = getAttenuationCoefficients(o.getMaterial())[e];
sum += att * lens[j];
}
// TODO: Normalization is never considered in the backprojectors,
// thus, iteratively applying forward and backward projections
// would yield to a scaling issue!
//
// length is in [mm]
// attenuation is in [g/cm^3]
// conversion from [g*mm/mc^3] = [g*0.1cm/cm^3] to [g/cm^2]
// --> sum/10.0;
double afterAttenuation = photonFlux[e] * Math.exp(-sum/10);
if (noise) {
double photonsWithNoise = StatisticsUtil.poissonRandomNumber(afterAttenuation);
if (energyIntegrating){
photonsWithNoise *= energies[e];
}
intensity += photonsWithNoise;
} else {
if (energyIntegrating){
afterAttenuation *= energies[e];
}
intensity += afterAttenuation;
}
}
return intensity;
}
/**
*
* @return the total intensity of the input spectrum
*/
public double getTotalIntensity() {
return inputSpectrum.getTotalIntensity();
}
/**
*
* @return the total photon flux of the input spectrum
*/
public double getTotalPhotonFlux() {
return inputSpectrum.getTotalPhotonFlux();
}
/**
*
* @return the time current product of the input spectrum
*/
public double getTimeCurrentProduct(){
return inputSpectrum.getTimeCurrentProduct();
}
/**
* @return the att
*/
public AttenuationType getAtt() {
return att;
}
/**
* @param att the att to set
*/
public void setAtt(AttenuationType att) {
this.att = att;
}
public double [] getAttenuationCoefficients(Material mat){
double [] result = attenuationCoefficientsMap.get(mat);
if (result == null){
MaterialUtils.updateAttenuationCoefficientsMap(attenuationCoefficientsMap, energies, mat, att);
result = attenuationCoefficientsMap.get(mat);
}
return result;
}
}