# coding: utf-8
# /*##########################################################################
#
# Copyright (c) 2016 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ###########################################################################*/
"""Create the link between freeart and a configuration file. Basically set up
a reconstruction
"""
__authors__ = ["H. Payno"]
__license__ = "MIT"
__date__ = "23/10/2017"
from freeart.configuration import config, structs
import freeart
import os
import sys
import freeart.utils
import freeart.utils.reconstrutils as reconstrutils
import freeart.utils.physicalelmts as physicalelmts
from freeart.utils import sinogramselection
import numpy as np
import logging
import fisx
from threading import Thread
from fisx import Detector
from fabio.edfimage import edfimage
from freeart.configuration import read
import numpy
_logger = logging.getLogger("freeARTInterpreter")
logging.basicConfig(level=logging.INFO)
elementInstance = None
[docs]class AbsConfInterpreter(object):
"""Abstract class from which each interpreter (fluorescence, tx...)
should inherit"""
def __init__(self, filePath, _config):
assert not (filePath is None and _config is None)
self.filePath = filePath
if filePath:
self.config = read(self.filePath)
else:
self.config = _config
assert(self.config.precision in config._ReconsConfig.PRECISION_TO_TYPE)
self.precisionType = config._ReconsConfig.PRECISION_TO_TYPE[self.config.precision]
self.setTheoreticalSampSize()
self.setLimitParallelReconstruction(-1)
self.setLimitVoxels(-1)
def setTheoreticalSampSize(self):
raise NotImplementedError('AbsConfInterpreter is a pure virtual class')
[docs] def iterate(self):
"""
Iterate on the reconstruction. Each iteration generate a backward
projection for each angle of acquisition (in a random order).
"""
raise NotImplementedError('AbsConfInterpreter is a pure virtual class')
def _reduceSinogram(self, matrix):
"""Reduce the data according to the DefinitionReductionfactor and the
ProjectionReductionFactor
:param matrix: the matrix to reduce
:return: the matrix reduced
"""
if self.config.definitionReduction is None:
self.config.definitionReduction = 1
matrix = self._reduceXData(matrix, self.config.definitionReduction)
if self.config.projectionNumberReduction is None:
self.config.projectionNumberReduction = 1
return self._reduceNbAnglesData(matrix, self.config.projectionNumberReduction)
def _reduceMatriceSize(self, data):
"""Reduce the matrice in x and y fron the setted configuration file
:param data: the matrix to reduce.
:return: a reduced matrix
"""
if data.shape[1] != self._getTheoreticalSampSize():
if self.config.definitionReduction is None:
self.dataXReduction = self.config.definitionReduction
res = None
reduction = self.config.definitionReduction
if len(data.shape) == 3:
res = data[:, ::reduction, ::reduction]
if len(data.shape) == 2:
res = data[::reduction, ::reduction]
if res is None:
raise RuntimeError("Unvalid matrix dimension")
else:
assert(res.shape[1] == self._getTheoreticalSampSize())
return res
else:
return data
def _getTheoreticalSampSize(self):
# return the size of the reconstructed sample
return self._thSampSize
def _getFirstFluoSino(self):
for sinoFile in self.config.sinograms:
for sino in self.config.sinograms[sinoFile]:
return sino
def _getAngles(self, nbAngles, selection):
"""
Compute nbAngle homogeneous angle between minAngle and maxAngle.
Take into account if the last angle should be equal to the first angle
from the config file (LastAngleEqualFirstAngle option)
:param nbAngles: the number of angle we want to generate
:return: the angle for homogeneous angle between each acquisition
"""
def getMinAngle():
if self.config.minAngle is None:
err = "can't find information about min angle in the .cfg file"
raise IOError(err)
return self.config.minAngle
def getMaxAngle():
if self.config.maxAngle is None:
err = "can't find information about max angle in the .cfg file"
raise IOError(err)
return self.config.maxAngle
# if specify that the last angle in the sinogram == first angle in the
# sinogram
if self.config.includeLastProjection is None:
self.config.includeLastProjection = False
minAngle = getMinAngle()
maxAngle = getMaxAngle()
if minAngle > maxAngle:
raise RuntimeError("minAngle > maxAngle, can't process")
if minAngle < -2.1*np.pi or minAngle > 2.1*np.pi or maxAngle < -2.1*np.pi or maxAngle > 2.1*np.pi:
err = """Angles should be defined in radians. We for minAngle and
maxAngle at least one angle < -2*pi or > 2.pi """
raise RuntimeError(err)
return sinogramselection.getAngles(minAngle=minAngle,
maxAngle=maxAngle,
selection=selection,
acquiInverted=self.config.acquiInverted,
lastAngleEqualFirst=self.config.includeLastProjection,
nbAngles=nbAngles)
def _reduceNbAnglesData(self, data, factorReduction):
"""
Reduce the number of angles in the data of a factor factorReduction
:param data: the matrix to reduce
:param factorReduction: the factor of reduction to apply
:return: the matrix reduced
.. warning:: the matrix should be a 2D numpy array
"""
return data[::factorReduction, ]
def _reduceXData(self, data, factorReduction):
"""
Reduce the number of projection in the data of a factor factorReduction
:param data: the matrix to reduce
:param factorReduction: the factor of reduction to apply
:return: the matrix reduced
.. warning:: the matrix should be a 2D numpy array
"""
return data[:, ::factorReduction]
def _centerSinogram(self, data, center):
"""Simple function to center the sinogram
:param data: the sinogram to center
:param center: The new center of rotation
"""
assert(len(data.shape) == 2)
width = min(abs(data.shape[1]-center), center)
if width > 0:
return data[:, center-width:center+width]
else:
_logger.info("can t center the sinogram. Centering avoid")
return data
[docs] def getSinograms(self):
"""
:return: A dictionnary of the sinograms to treat"""
err = 'getSinograms not implemented. '
err += 'AbsConfInterpreter is a pure virtual class'
raise NotImplementedError(err)
[docs] def getFile(self):
"""
:return: the current configuratioin file """
return self.filePath
[docs] def setLimitParallelReconstruction(self, val):
"""
Limit the number of parallel reconstruction we can run
:param val: the new limitation of parallel reconstructions to run
"""
if val > 0:
self.limitParaRecon = val
if val == -1:
self.limitParaRecon = sys.maxsize
[docs] def setLimitVoxels(self, val):
"""
Limit the maximal number of voxel we can run during one iteration
:param val: the maximal number of voxel to run during one step.
"""
if val > 0:
self.limitVoxel = val
if val == -1:
self.limitVoxel = sys.maxsize
def _normalizeRawData(self, sinogram, reduceI0Data=False):
"""Normalization of raw data. Apply the new center of rotation, reduce
data if needed and normalize fron I0 value.
:param sinogram: the sinogram to normalize
:param reduceI0Data: if true then normalize the sinogram size if
needed. Option needed when part of the data has
already been normalized (for example when a
transmission reconstruction has been run before).
"""
if self.config.useAFileForI0 is True and type(self.config.I0) is str:
sinoI0 = reconstrutils.LoadEdf_2D(self.I0)
# warning : we should always first center the sinogram then reduce
# is dimension
if self.config.center:
sinoI0 = self._centerSinogram(sinoI0, self.config.center)
if reduceI0Data:
sinoI0 = self._reduceSinogram(sinoI0)
if (sinoI0.shape != sinogram.shape):
err = """I0 sinogram (%s) and sinogram (%s) have different
shapes""" % (sinoI0.shape, sinogram.shape)
raise ValueError(err)
else:
sinogram = sinogram / sinoI0
elif self.config.I0 is not None:
if self.config.I0 <= 0.0:
raise ValueError("invalid I0 value given")
else:
sinogram = sinogram/self.config.I0
return sinogram
[docs]class TxConfInterpreter(AbsConfInterpreter):
"""
Interpreter for tx reconstruction
:param filePath: the path of the cfg file
"""
txRecons = None
def __init__(self, filePath, _config=None):
AbsConfInterpreter.__init__(self, filePath=filePath, _config=_config)
if _config:
assert(isinstance(_config, config.TxConfig))
# make sure data is here
self._setUpTxReconstruction()
[docs] def getReconstructionAlgorithms(self):
"""
:return: the dictionay of the freeart reconstruction.
"""
return {"TxReconstruction": self.txRecons}
def setTheoreticalSampSize(self):
assert(self.config.sinogram)
self.config.sinogram.loadData(refFile=self.filePath)
assert(self.config.sinogram.data is not None)
sino = self.config.sinogram
assert isinstance(sino, structs.TxSinogram)
mod = sino.data.shape[1] % self.config.definitionReduction
mod = 1 if mod > 0 else 0
self._thSampSize = sino.data.shape[1] / self.config.definitionReduction + mod
def _setUpTxReconstruction(self):
"""
Setup the FreeART Tx algorithm to fit with the self.configuration file
"""
self.txRecons = self._buildAlgoTx()
def _buildAlgoTx(self):
"""
Build a transmission ART algorithm for the given algorithm and given
self.configuration file
"""
# step 1 : normalization of the sinogram
if self.config.center:
self.config.sinogram.data = self._centerSinogram(self.config.sinogram.data, self.config.center)
# step 2 : deal with data reduction.
# Rule : do not reduce data when the target is to produce an absMat.
# Simplify a lot stuff (specially when loading fluo config file)
self.config.sinogram.data = self._reduceSinogram(self.config.sinogram.data)
# step 3 : get angles
selection = self.config.projections
self.angles = self._getAngles(nbAngles=self.config.sinogram.data.shape[0],
selection=selection)
# step 4 : normalize from I0
self.config.sinogram.data = self._normalizeRawData(self.config.sinogram.data,
reduceI0Data=True)
# reduce the sinogram if we want to avoid some projection
self.config.sinogram.data = sinogramselection.getSelection(
selection=selection,
projections=self.config.sinogram.data)
# step 5 : move it to 3D
self.config.sinogram.data.shape = (1,
self.config.sinogram.data.shape[0],
self.config.sinogram.data.shape[1])
# if some values of the sinogram are negative and the sinogram is I
# (and not -log(I))
if self.config.computeLog and (True in (self.config.sinogram.data < 0.0)):
err = """given sinogram has negative values. Can't process
the reconstruction"""
raise RuntimeError(err)
# create the adapted ARTAlgorithm
if self.config.computeLog is True:
logVal = -np.log(self.config.sinogram.data)
else:
logVal = self.config.sinogram.data
assert(np.nan not in logVal)
txAlgo = freeart.TxBckProjection(logVal.astype(self.precisionType),
self.angles)
# step 6 : setting up extra parameters
if self.config.voxelSize:
txAlgo.setVoxelSize(self.config.voxelSize * self.config.definitionReduction)
if self.config.oversampling:
txAlgo.setOverSampling(self.config.oversampling)
if self.config.I0 and not (self.config.useAFileForI0 is True and self.config.I0):
txAlgo.setI0(self.config.I0)
if self.config.beamCalcMethod:
txAlgo.setRayPointCalculationMethod(self.config.beamCalcMethod)
if self.config.dampingFactor:
txAlgo.setDampingFactor(self.config.dampingFactor)
return txAlgo
[docs] def getSinograms(self):
"""
:return: A dictionnary of the sinograms to treat
"""
return [reconstrutils.decreaseMatSize(self.config.sinogram.data.copy())]
[docs] def iterate(self, nbIteration):
"""
Iterate on the reconstruction. Each iteration generate a backward
projection for each angle of acquisition (in a random order).
"""
if self.txRecons is None:
err = "No ART algorithm defined"
raise ValueError(err)
return self.txRecons.iterate(nbIteration)
def setRandSeedToZero(self, b):
if self.txRecons is None:
err = "No ART algorithm defined"
raise ValueError(err)
self.txRecons.setRandSeedToZero(b)
[docs]class FluoConfInterpreter(AbsConfInterpreter):
"""
FreeART interpreter for fluorescence reconstruction.
Basicall y take a .cfg file containing the description of the
reconstruction in input and create all corresponding ARTAlgorithm for such
reconstructions
:param filePath: the configuration file to be interpreted.
"""
def __init__(self, filePath, _config=None):
AbsConfInterpreter.__init__(self, filePath=filePath, _config=_config)
if _config:
assert (isinstance(_config, config.FluoConfig))
self.materials = {}
"""List all the materials present in the sample"""
self.selfAbsMatFile = {}
"""File where to get the selfAbsMat file in case the user give them"""
self.absMatrixFile = None
"""The file of the absorption matrix file"""
self.fluoSinogramsAlgorithms = {}
"all the art algorithm to fit the configuration file"
self.fluoSinograms = {}
"sinograms to treat"
self.interactionMatrice = {}
"the interaction matrices"
self.angles = None
self._sinoElemt = {}
self._setUpFluoReconstruction()
self._buildAlgorithms()
def setTheoreticalSampSize(self):
sino = self._getFirstFluoSino()
assert(sino)
assert isinstance(sino, structs.FluoSino)
mod = sino.data.shape[1] % self.config.definitionReduction
mod = 1 if mod > 0 else 0
self._thSampSize = sino.data.shape[1] / self.config.definitionReduction + mod
def _setUpFluoReconstruction(self):
"""
Setup the FreeART Fluo algorithm to fit with the self.configuration
file
"""
if len(self.config.sinograms) < 1:
raise ValueError("No fluorescence sinograms given")
if self.config.absMat is None or (self.config.absMat.fileInfo is None and self.config.absMat.data is None):
raise ValueError('Absorption matrix not defined')
if self.config.absMat.data is None and not os.path.isfile(self.config.absMat.fileInfo.getFile()):
err = self.config.absMat + " absorption file setted does not exist"
raise ValueError(err)
# in the case we need to build the absorption matrice from ourself
if self.config.isAbsMatASinogram is True:
mess = """The interpreter is not able to reconstruct by himself the
absorption matrix. Please build it first, then give if back to
the freeart interpreter"""
raise RuntimeError(mess)
elif self.config.absMat.data is None:
self.config.absMat.loadData(self.filePath)
self.config.absMat.data = self._reduceMatriceSize(self.config.absMat.data)
self.config.absMat.data.shape = (self.config.absMat.data.shape[0],
self.config.absMat.data.shape[1],
1)
self.__prepareSelfAbs()
def __prepareMaterials(self):
"""
prepare materials for fisx
"""
if self.config.materials is None or self.config.materials.materials is None or \
self.config.materials.matComposition is None or \
self.config.materials.materials.fileInfo is None or \
self.config.materials.matComposition.fileInfo is None:
self.materials = None
return
else:
assert(self.config.materials.materials is not None)
assert(self.config.materials.materials.fileInfo is not None)
assert(self.config.materials.matComposition is not None)
assert(self.config.materials.matComposition.fileInfo is not None)
assert(os.path.isfile(self.config.materials.materials.fileInfo.filePath))
self.config.materials.loadData(self.filePath)
materials = self.config.materials.materials.data
physicalelmts.elementInstance.removeMaterials()
for materialName in materials:
self.materials[materialName] = reconstrutils.convertToFisxMaterial(
materialName,
materials[materialName])
physicalelmts.elementInstance.addMaterial(
self.materials[materialName])
def __getFisxDetector(self):
"""
:return: the fisx detector to compute the interaction matrix
"""
# TODO : get the material of the detector
# pb : in this case we will take twice the effect of the solid angle.
# Something that we don't want to...
detector = Detector("Si1", 2.33, 1.0)
detector.setDiameter(self.config.detector.width)
@staticmethod
[docs] def getProbabilityOfEmission(fisxMat, energy, element, fisxElements,
shell=None):
"""
compute the probability of emission of the element for a specific
energy and material.
:param fisxMat: the fisx material of the interaction
:param energy: the nergy of the incomng beam
:param element: the generated element
:param shell: the targetted shell of the interaction.
:return: the probability of emission of the element for a specific
energy and material.
"""
if not type(fisxMat) is fisx.Material:
raise TypeError('material should be a fisx.Material')
excitations = fisxElements.getExcitationFactors(element, [energy])
if sys.version_info[1] > 2:
photonelectricMassAttenuation = fisxElements.getMassAttenuationCoefficients(fisxMat.getComposition(), energy)[b'photoelectric']
else:
photonelectricMassAttenuation = fisxElements.getMassAttenuationCoefficients(fisxMat.getComposition(), energy)['photoelectric']
if shell is None:
rate = 0.
for excitation in excitations:
rate += excitations[excitation]['rate'] * photonelectricMassAttenuation[0]
return rate
else:
if shell not in excitations:
raise ValueError('Required shell is not existing for the elemnt')
else:
return excitations[element]['rate'][shell] * photonelectricMassAttenuation
def __prepareSelfAbs(self):
"""
Function to prepare the self absorption matrix
"""
self.__prepareMaterials()
self.selfAbsMatrixNeedToBEComputed = self.materials is not None
if self.selfAbsMatrixNeedToBEComputed is True:
# fluo case
self.E0 = self.config.E0
self.sampleComposition = self.config.materials.matComposition.data
def _resetFluoSinograms(self):
"""Reset all fluorescence sinograms """
self.fluoSinogramsAlgorithms = {}
self.fluoSinograms = {}
self.interactionMatrice = {}
def _resetAngles(self):
"""Reset all angles """
self.angles = None
def _resetSinoElemt(self):
"""
Reset the sinogram element
"""
self._sinoElemt = {}
def _buildAlgorithms(self):
"""
will setup one ART algorithm per fluo sinogram
"""
self._resetFluoSinograms()
self._resetAngles()
[self._buildAlgorithm(filePath) for filePath in self.config.sinograms]
def __getTitle(self, header):
"""return the title of the header"""
if 'title' in header:
return header['title']
if 'Title' in header:
return header['Title']
return None
def __getMeanData(self, edfReader, groupFiles, indexDataToTreat, referenceName):
"""Simply create a mean sinogram from a list of files and check that
each sinogram of the file as the same ID (referenceName)
:param edfReader: the edfReader for the file
:param groupFiles: the list of file definig the group
:param indexDataToTreat: index of the sinogram in the file
:param referenceName: The 'title' of the reference sinogram
"""
sino = None
for file in groupFiles:
# check that titles are the same as the reference file. otherwise
# through an error
if edfReader.nframes == 1:
if self.__getTitle(edfReader.getHeader()) != referenceName:
raise ValueError('incoherent title name between sinogram to compute the mean')
if sino is None:
sino = edfReader.getData()
else:
sino += edfReader.getData()
else:
if self.__getTitle(edfReader.getframe(int(indexDataToTreat)).getHeader()) != referenceName:
raise ValueError('incoherent title name between sinogram to compute the mean')
if sino is None:
sino = edfReader.getframe(int(indexDataToTreat)).getData()
else:
sino += edfReader.getframe(int(indexDataToTreat)).getData()
return sino / len(groupFiles)
def _buildAlgorithm(self, filePath):
"""build the algorithm for the givem file path"""
# read the file
if filePath is None:
err = "unvalid file " + self.filePath
raise IOError(err)
if not os.path.isfile(filePath):
err = "file " + filePath + " not existing in your system"
raise ValueError(err)
# extract sinograms
assert(filePath in self.config.sinograms)
datasets = self.config.sinograms[filePath]
gf, groupType = reconstrutils.tryToFindPattern(filePath)
if groupType == 'pymca':
try:
self.__treatAsAPymcaGroupFiles(filePath, datasets)
except ValueError:
err = """Found incoherence in the files or in sinograms, unable
to associate pymca group files"""
raise ValueError(err)
else:
self.__treatAsAStandardFile(filePath, datasets)
def __treatAsAPymcaGroupFiles(self, filePath, sinograms):
"""
create all the algorithm from getting the mean file from the pymca
"""
groupFiles = reconstrutils.tryToFindPattern(filePath, 'pymca')
assert(len(groupFiles) > 0)
if (len(groupFiles) < 2):
err = "file can't be associated with any other pymca file"
raise ValueError(err)
edfReader = edfimage()
edfReader.read(filePath)
for sino in sinograms:
assert(isinstance(sino, structs.FluoSino))
name = sino.name
self._sinoElemt[name] = sino.physElmt
EF = sino.EF
refName = None
if edfReader.nframes == 1:
refName = self.__getTitle(edfReader.getHeader())
else:
refName = self.__getTitle(
edfReader.getframe(int(sino.index)).getHeader())
sino.data = self.__getMeanData(edfReader=edfReader,
groupFiles=groupFiles,
indexDataToTreat=sino.index,
referenceName=refName)
self._addFluoARTAlgorithm(sino)
def __treatAsAStandardFile(self, filePath, datasets):
"""
create all the algorithm at the standard way
"""
edfReader = edfimage()
edfReader.read(filePath)
for filePath in self.config.sinograms:
for sino in self.config.sinograms[filePath]:
assert(isinstance(sino, structs.FluoSino))
self._sinoElemt[sino.name] = sino.physElmt
sino.loadData(refFile=filePath)
self._addFluoARTAlgorithm(sino)
def _updateSelfAbsMat(self, fluoSinogram):
"""
.. note:: __prepareSelfAbs need to be called before
"""
assert(isinstance(fluoSinogram, structs.FluoSino))
if self.config.reconsType == 'Compton':
# compton case
return self.absMatrix
elif self.selfAbsMatrixNeedToBEComputed is False:
if fluoSinogram.selfAbsMat is None:
raise ValueError('Can\'t build the required algorithm, selfAbsMat should be defined')
fluoSinogram.selfAbsMat.loadData(self.filePath)
return fluoSinogram.selfAbsMat.data
else:
assert(self.config.reconsType == 'Fluorescence')
assert (isinstance(fluoSinogram, structs.FluoSino))
if fluoSinogram.EF is None:
err = "No EF given. Can't process the fluorescence reconstruction"
raise ValueError(err)
assert(fluoSinogram.EF > 0.0)
# In the case of the fluorescence
selfAbsMatrix = np.zeros(self.config.absMat.data.shape, dtype=np.float64)
for materialName in self.materials:
cross_section_EF_mat = physicalelmts.elementInstance.getElementMassAttenuationCoefficients(materialName, fluoSinogram.EF)
cross_section_E0_mat = physicalelmts.elementInstance.getElementMassAttenuationCoefficients(materialName, self.E0)
# Note : energy is in keV in freeart AND in physx. So no need for conversion
condition = np.in1d(self.sampleComposition.ravel(), materialName).reshape(selfAbsMatrix.shape)
selfAbsMatrix[condition] = self.config.absMat.data[condition] * cross_section_EF_mat['total'][0] / cross_section_E0_mat['total'][0]
# reduce the sample composition if needed (if sinogram definition
# is reduced )
assert(selfAbsMatrix.ndim == 3)
selfAbsMatrix = self._reduceMatriceSize(selfAbsMatrix)
assert(self.config.absMat.data.shape == selfAbsMatrix.shape)
return selfAbsMatrix
def _addFluoARTAlgorithm(self, fluoSinogram):
"""
Add an ART algorithm to the set of algorithm
:param fluoSinogram: the sinogram for which we want to build an ART
algorithm
"""
assert(isinstance(fluoSinogram, structs.FluoSino))
assert(fluoSinogram.data is not None)
_logger.info('create the algorithm for ' + fluoSinogram.physElmt)
if self.config.center:
fluoSinogram.data = self._centerSinogram(data=fluoSinogram.data,
center=self.config.center)
fluoSinogram.data = self._normalizeRawData(fluoSinogram.data)
fluoSinogram.data = self._reduceSinogram(fluoSinogram.data)
selection = self.config.projections
self.angles = self._getAngles(nbAngles=fluoSinogram.data.shape[0],
selection=selection)
# reduce to some projection
fluoSinogram.data = sinogramselection.getSelection(selection=selection,
projections=fluoSinogram.data)
# Do not do this, already managed in the freeart core process
# if self.angles[0] != 0.0:
# # rotate the detector position from the origon of a certain angle
# self.detPos = tuple(freeart.utils.computeDetectorPosition(
# detPos=self.detPos,
# angle=self.angles[0]))
# self.detSetup = [(self.detPos, self.detWidth)]
#
# info = 'Updating detector position for the first projection.'
# info += 'New position is %s' % str(self.detPos)
# _logger.info(info)
# move it to 3D
fluoSinogram.data.shape = (1, fluoSinogram.data.shape[0], fluoSinogram.data.shape[1])
# for now in compton Mode
if not fluoSinogram.data.shape[2] == self.config.absMat.data.shape[1]:
error = "can't run the reconstruction, width of the sinogram (%s)" \
"and the absorption matrix dim (%s) are differente" \
% (fluoSinogram.data.shape[2], self.config.absMat.data.shape[1])
raise ValueError(error)
selfAbsMatrix = self._updateSelfAbsMat(fluoSinogram)
self.config.absMat.data = self.__checkAbsMatDim(self.config.absMat.data)
selfAbsMatrix = self.__checkAbsMatDim(selfAbsMatrix)
if not self.config.absMat.data.shape == selfAbsMatrix.shape:
error = "can't run the reconstruction, shapes of the absorption " \
"matrix (%s )" % str(self.config.absMat.data.shape)
error += "and shape of the self absorption matrix are not fitting" \
" (%s )" % str(selfAbsMatrix.shape)
_logger.error(error)
raise ValueError(error)
# creating the ART algorithm
detector = self.config.detector
detSetup = [(detector.getPosition(), detector.width)]
alRecons = freeart.FluoBckProjection(sinoDat=fluoSinogram.data.astype(self.precisionType),
sinoAngles=self.angles,
expSetUp=detSetup,
absorp=self.config.absMat.data.astype(self.precisionType),
selfAbsorp=selfAbsMatrix.astype(self.precisionType))
if self.config.voxelSize:
alRecons.setVoxelSize(self.config.voxelSize * self.config.definitionReduction)
if self.config.oversampling:
alRecons.setOverSampling(self.config.oversampling)
if self.config.I0:
alRecons.setI0(self.config.I0)
if self.config.beamCalcMethod:
alRecons.setRayPointCalculationMethod(self.config.beamCalcMethod)
if self.config.outBeamCalMethod:
alRecons.setOutgoingRayAlgorithm(self.config.outBeamCalMethod)
if self.config.solidAngleOff:
alRecons.turnOffSolidAngle(self.config.solidAngleOff)
if self.config.dampingFactor:
alRecons.setDampingFactor(self.config.dampingFactor)
self.fluoSinogramsAlgorithms[fluoSinogram.physElmt] = alRecons
self.fluoSinograms[fluoSinogram.physElmt] = fluoSinogram
# create the interaction matrice
self.interactionMatrice[fluoSinogram.physElmt] = self.__computeInteractionMatrice(fluoSinogram.physElmt)
def __checkAbsMatDim(self, mat):
if mat.ndim != 3:
assert(mat.ndim is 2)
mat = mat.reshape(mat.shape[0], mat.shape[1], 1)
return mat
def __computeInteractionMatrice(self, name):
"""Return the interaction matrice according to the materials in the
sample and the element
:param name: the id of the element for which we want to compute the
interaction matrix
"""
if self.config.reconsType == 'Compton':
interactionMatrice = np.ones((self.config.absMat.data.shape[0],
self.config.absMat.data.shape[1]),
dtype=np.float64)
else:
if hasattr(self, 'materials'):
interactionMatrice = np.zeros((self.config.absMat.data.shape[0],
self.config.absMat.data.shape[1]),
dtype=np.float64)
if self.materials is not None:
for materialName in self.materials:
probabilityOfEmission = FluoConfInterpreter.getProbabilityOfEmission(
fisxMat=self.materials[materialName],
energy=self.config.E0,
element=self._sinoElemt[name],
fisxElements=physicalelmts.elementInstance,
shell=None)
condition = np.in1d(self.sampleComposition.ravel(),
materialName).reshape(interactionMatrice.shape)
interactionMatrice[condition] = probabilityOfEmission
interactionMatrice.shape = (interactionMatrice.shape[0],
interactionMatrice.shape[1],
1)
else:
mes = """no materials specified, set an interaction matrice to
ones"""
_logger.warning(mes)
interactionMatrice = np.ones((self.config.absMat.data.shape[0],
self.config.absMat.data.shape[1]),
dtype=np.float64)
return interactionMatrice
[docs] def getReconstructionAlgorithms(self):
"""
:return: all the algorithm created for reconstructions
"""
return self.fluoSinogramsAlgorithms
[docs] def getSinograms(self):
"""
:return: a dictionnary with all the sinograms we are trying to
reconstruct
"""
return self.fluoSinograms
[docs] def iterate(self, nbIteration):
"""
Iterate on the reconstruction. Each iteration generate a backward
projection for each angle of acquisition (in a random order).
"""
threads = []
nbVoxel = 0
for index, (algoID, algo) in enumerate(self.fluoSinogramsAlgorithms.items()):
if index % self.limitParaRecon == 0 or (nbVoxel > self.limitVoxel):
[thread.start() for thread in threads]
[thread.join() for thread in threads]
threads = []
assert(not self.fluoSinogramsAlgorithms[algoID] is None)
threads.append(Thread(target=self.fluoSinogramsAlgorithms[algoID].iterate,
args=(nbIteration,)))
nbVoxel = self.fluoSinograms[algoID].data.shape[1] * self.fluoSinograms[algoID].data.shape[1]
[thread.start() for thread in threads]
[thread.join() for thread in threads]
threads = []
def setRandSeedToZero(self, b):
for elmt in self.fluoSinogramsAlgorithms:
assert(self.fluoSinogramsAlgorithms[elmt] is not None)
self.fluoSinogramsAlgorithms[elmt].setRandSeedToZero(b)
[docs] def getDensityPhantoms(self):
"""
Return the phantom taking into account the probability of generation
of the physical element
"""
res = {}
for element in self.fluoSinogramsAlgorithms:
ph = np.zeros(self.interactionMatrice[element].shape)
if self.interactionMatrice[element] is not None:
probInter = self.fluoSinogramsAlgorithms[element].getPhantom()[self.interactionMatrice[element] > 0.0] / \
self.interactionMatrice[element][self.interactionMatrice[element] > 0.0]
ph[self.interactionMatrice[element] > 0.0] = probInter
res[element] = ph
return res
def getInteractionPhantoms(self):
return self.interactionMatrice