--- /dev/null
+#=========================================================
+# Install scripts when running make install
+SET(SCRIPTS
+ gate_job_cluster.job
+ gate_make_merge_release.sh
+ gate_make_release.sh
+ gate_power_merge.sh
+ gate_run_submit_cluster.sh
+ gate_upload_release.sh
+ mergeDosePerEnergyFile.sh
+ mergeStatFile.py
+ mergeStatFile.sh
+)
+
+INSTALL (FILES ${SCRIPTS} DESTINATION ${CMAKE_INSTALL_PREFIX}/bin PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_EXECUTE WORLD_EXECUTE)
--- /dev/null
+#!/bin/bash
+
+set -u
+set -e
+
+function error {
+echo "ERROR: $1"
+echo "$(basename $0)"
+exit 1
+}
+
+function get_deps {
+targetbinary=${1}
+targetdir=${2}
+test -d ${targetdir} || error "${targetdir} isn't a directory"
+ldd ${targetbinary} | while read library; do
+ libfile="$(echo ${library} | awk -F' ' '/=> \// {print $3}')"
+ test $libfile || continue # didn't macht regex
+ test -f "${targetdir}/$(basename ${libfile})" && continue # already exists
+ cp "${libfile}" "${targetdir}"
+ get_deps "${libfile}" "${targetdir}"
+done
+}
+
+function package_target {
+targetname="${1}"
+which "${targetname}" > /dev/null || error "cant locate ${targetname}"
+targetbin="$(which ${targetname})"
+echo "${targetname} executable is ${targetbin}"
+
+echo "Getting libraries"
+targetdir="$(mktemp -d)"
+get_deps "${targetbin}" "${targetdir}"
+
+echo "Removing unused libraries"
+rm -f ${targetdir}/libdl.so*
+rm -f ${targetdir}/libm.so*
+rm -f ${targetdir}/libstdc++.so*
+rm -f ${targetdir}/libgcc_s.so*
+rm -f ${targetdir}/libpthread.so*
+rm -f ${targetdir}/libc.so*
+
+echo "Copying binary"
+cp "${targetbin}" .
+for filename in $(find ${targetdir} -name '*.so*'); do cp ${filename} . ; done
+
+echo "Cleaning up"
+rm -r "${targetdir}"
+}
+
+filenames=("clitkImageArithm" "clitkMergeRootFiles" "clitkMergeAsciiDoseActor" "clitkImageUncertainty" "mergeStatFile.sh" "gate_power_merge.sh")
+
+for input in "${filenames[@]}"; do
+ package_target "${input}" || error "error while packaging ${input}"
+done
+
+echo "Making release"
+tar -czvf merge_release.tar.gz ** \
+ || usage "can't create release zip"
+
+
+
#!/usr/bin/env bash
-set -u
+set -u
function error {
-echo "ERROR: $1"
-exit 1
+ echo "ERROR: $1"
+ exit 1
}
warning_count=0
function warning {
-let "warning_count++"
-echo "MERGE_WARNING: $1"
+ let "warning_count++"
+ echo "MERGE_WARNING: $1"
}
function start_bar {
-count_max="${1:?"provide count max"}"
+ count_max="${1:?"provide count max"}"
}
function update_bar {
-local count="${1:?"provide count"}"
-local message="${2:?"provide message"}"
-local percent=$(echo "100*${count}/${count_max}" | bc)
-#printf "[%03d/%03d] %3d%% %-80.80s\r" ${count} ${count_max} ${percent} "${message}"
-printf "[%03d/%03d] %3d%% %-80.80s\n" ${count} ${count_max} ${percent} "${message}"
+ local count="${1:?"provide count"}"
+ local message="${2:?"provide message"}"
+ local percent=$(echo "100*${count}/${count_max}" | bc)
+ #printf "[%03d/%03d] %3d%% %-80.80s\r" ${count} ${count_max} ${percent} "${message}"
+ printf "[%03d/%03d] %3d%% %-80.80s\n" ${count} ${count_max} ${percent} "${message}"
}
function end_bar {
-unset count_max
-#echo -ne '\n'
+ unset count_max
+ #echo -ne '\n'
}
function check_interfile {
-local input_interfile="${1:?"provide input interfile"}"
+ local input_interfile="${1:?"provide input interfile"}"
-grep -qs '!INTERFILE :=' "${input_interfile}" || return 1
+ grep -qs '!INTERFILE :=' "${input_interfile}" || return 1
-local header_byte_size=$(awk -F' ' '
+ local header_byte_size=$(awk -F' ' '
BEGIN { zsize = 0; }
/matrix size/ && $3 == "[1]" { xsize = $5; }
/matrix size/ && $3 == "[2]" { ysize = $5; }
/number of bytes per pixel/ { byte_per_pixel = $7; }
END { print xsize * ysize * zsize * byte_per_pixel; }' "${input_interfile}")
-local raw_interfile="$(dirname "${input_interfile}")/$(awk -F' := ' '/name of data file/ { print $2; }' "${input_interfile}")"
+ local raw_interfile="$(dirname "${input_interfile}")/$(awk -F' := ' '/name of data file/ { print $2; }' "${input_interfile}")"
-test -f "${raw_interfile}" || return 1
-test $(stat -c%s "${raw_interfile}") -eq ${header_byte_size} || return 1
+ test -f "${raw_interfile}" || return 1
+ test $(stat -c%s "${raw_interfile}") -eq ${header_byte_size} || return 1
}
function write_mhd_header {
-local input_interfile="${1:?"provide input interfile"}"
-local output_mhd="$(dirname "${input_interfile}")/$(basename "${input_interfile}" ".hdr").mhd"
+ local input_interfile="${1:?"provide input interfile"}"
+ local output_mhd="$(dirname "${input_interfile}")/$(basename "${input_interfile}" ".hdr").mhd"
-check_interfile "${input_interfile}" || error "${input_interfile} isn't an interfile image"
+ check_interfile "${input_interfile}" || error "${input_interfile} isn't an interfile image"
-local header_start='ObjectType = Image
+ local header_start='ObjectType = Image
NDims = 3
AcquisitionDate = none
BinaryData = True
DistanceUnits = mm
AnatomicalOrientation = RIP'
-echo "${header_start}" > "${output_mhd}"
+ echo "${header_start}" > "${output_mhd}"
-awk -F' ' '
+ awk -F' ' '
/scaling factor/ && $4 == "[1]" { xspacing = $6; }
/scaling factor/ && $4 == "[2]" { yspacing = $6; }
END { print "ElementSpacing = " xspacing " " yspacing " 1"; }' "${input_interfile}" >> "${output_mhd}"
-awk -F' ' '
+ awk -F' ' '
BEGIN { zsize = 0; }
/matrix size/ && $3 == "[1]" { xsize = $5; }
/matrix size/ && $3 == "[2]" { ysize = $5; }
/number of projections/ { zsize += $5; }
END { print "DimSize = " xsize " " ysize " " zsize; }' "${input_interfile}" >> "${output_mhd}"
-awk -F' := ' '
+ awk -F' := ' '
/number format/ { format = $2; }
/number of bytes per pixel/ { byte_per_pixel = $2 }
END {
print "ElementType = MET_INT";
}' "${input_interfile}" >> "${output_mhd}"
-awk -F' := ' '
+ awk -F' := ' '
/name of data file/ { print "ElementDataFile = " $2; }' "${input_interfile}" >> "${output_mhd}"
}
test -x "./clitkMergeRootFiles" && rootMerger="./clitkMergeRootFiles"
function merge_root {
-local merged="$1"
-shift
-echo " ${indent}entering root merger"
-echo " ${indent}merger is ${rootMerger}"
-echo " ${indent}creating ${merged}"
-#echo "######## $#"
-#echo "######## $*"
-
-if test $# -eq 1
-then
- echo " ${indent}just one partial file => just copy it"
- cp "$1" "${merged}"
- return
-fi
-
-local count=0
-local arguments=" -o ${merged}"
-while test $# -gt 0
-do
- local partial="$1"
+ local merged="$1"
shift
- let count++
- local arguments=" -i ${partial} ${arguments}"
-done
-${rootMerger} ${arguments} > /dev/null || error "error while calling ${rootMerger}"
-echo " ${indent}merged ${count} files"
+ echo " ${indent}entering root merger"
+ echo " ${indent}merger is ${rootMerger}"
+ echo " ${indent}creating ${merged}"
+ #echo "######## $#"
+ #echo "######## $*"
+
+ if test $# -eq 1
+ then
+ echo " ${indent}just one partial file => just copy it"
+ cp "$1" "${merged}"
+ return
+ fi
+
+ local count=0
+ local arguments=" -o ${merged}"
+ while test $# -gt 0
+ do
+ local partial="$1"
+ shift
+ let count++
+ local arguments=" -i ${partial} ${arguments}"
+ done
+ ${rootMerger} ${arguments} > /dev/null || error "error while calling ${rootMerger}"
+ echo " ${indent}merged ${count} files"
}
statMerger="mergeStatFile.py"
test -x "./mergeStatFile.sh" && statMerger="./mergeStatFile.sh"
function merge_stat {
-local merged="$1"
-shift
-echo " ${indent}entering stat merger"
-echo " ${indent}merger is ${statMerger}"
-echo " ${indent}creating ${merged}"
-local count=0
-start_bar $#
-while test $# -gt 0
-do
- local partial="$1"
+ local merged="$1"
shift
- let count++
-
- if test ! -f "${merged}"
- then
- update_bar ${count} "copying first partial result ${partial}"
- cp "${partial}" "${merged}"
- continue
- fi
+ echo " ${indent}entering stat merger"
+ echo " ${indent}merger is ${statMerger}"
+ echo " ${indent}creating ${merged}"
+ local count=0
+ start_bar $#
+ while test $# -gt 0
+ do
+ local partial="$1"
+ shift
+ let count++
+
+ if test ! -f "${merged}"
+ then
+ update_bar ${count} "copying first partial result ${partial}"
+ cp "${partial}" "${merged}"
+ continue
+ fi
- update_bar ${count} "adding ${partial}"
- ${statMerger} -i "${merged}" -j "${partial}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${statMerger}"
-done
-end_bar
-echo " ${indent}merged ${count} files"
+ update_bar ${count} "adding ${partial}"
+ ${statMerger} -i "${merged}" -j "${partial}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${statMerger}"
+ done
+ end_bar
+ echo " ${indent}merged ${count} files"
}
doseMerger="mergeDosePerEnegryFile.sh"
test -x "./mergeDosePerEnergyFile.sh" && doseMerger="./mergeDosePerEnergyFile.sh"
function merge_dose {
-local merged="$1"
-shift
-echo " ${indent}entering dose merger"
-echo " ${indent}merger is ${doseMerger}"
-echo " ${indent}creating ${merged}"
-local count=0
-start_bar $#
-while test $# -gt 0
-do
- local partial="$1"
+ local merged="$1"
shift
- let count++
-
- if test ! -f "${merged}"
- then
- update_bar ${count} "copying first partial result ${partial}"
- cp "${partial}" "${merged}"
- continue
- fi
+ echo " ${indent}entering dose merger"
+ echo " ${indent}merger is ${doseMerger}"
+ echo " ${indent}creating ${merged}"
+ local count=0
+ start_bar $#
+ while test $# -gt 0
+ do
+ local partial="$1"
+ shift
+ let count++
+
+ if test ! -f "${merged}"
+ then
+ update_bar ${count} "copying first partial result ${partial}"
+ cp "${partial}" "${merged}"
+ continue
+ fi
- update_bar ${count} "adding ${partial}"
- ${doseMerger} -i "${merged}" -j "${partial}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${doseMerger}"
-done
-end_bar
-echo " ${indent}merged ${count} files"
+ update_bar ${count} "adding ${partial}"
+ ${doseMerger} -i "${merged}" -j "${partial}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${doseMerger}"
+ done
+ end_bar
+ echo " ${indent}merged ${count} files"
}
txtImageMerger="clitkMergeAsciiDoseActor"
test -f "./clitkMergeAsciiDoseActor" && txtImageMerger="./clitkMergeAsciiDoseActor"
function merge_txt_image {
-local merged="$1"
-shift
-echo " ${indent}entering text image merger"
-echo " ${indent}merger is ${txtImageMerger}"
-echo " ${indent}creating ${merged}"
-local count=0
-start_bar $#
-while test $# -gt 0
-do
- local partial="$1"
+ local merged="$1"
shift
- let count++
-
- if test ! -f "${merged}"
- then
- update_bar ${count} "copying first partial result ${partial}"
- cp "${partial}" "${merged}"
- continue
- fi
+ echo " ${indent}entering text image merger"
+ echo " ${indent}merger is ${txtImageMerger}"
+ echo " ${indent}creating ${merged}"
+ local count=0
+ start_bar $#
+ while test $# -gt 0
+ do
+ local partial="$1"
+ shift
+ let count++
+
+ if test ! -f "${merged}"
+ then
+ update_bar ${count} "copying first partial result ${partial}"
+ cp "${partial}" "${merged}"
+ continue
+ fi
- update_bar ${count} "adding ${partial}"
- local header="$(cat "${merged}" | head -n 6)"
- local tmp="$(mktemp)"
- ${txtImageMerger} -i "${partial}" -j "${merged}" -o "${tmp}" 2> /dev/null > /dev/null || warning "error while calling ${txtImageMerger}"
- echo "${header}" > "${merged}"
- grep -v '## Merge' "${tmp}" >> "${merged}"
- rm "${tmp}"
-done
-end_bar
-echo " ${indent}merged ${count} files"
+ update_bar ${count} "adding ${partial}"
+ local header="$(cat "${merged}" | head -n 6)"
+ local tmp="$(mktemp)"
+ ${txtImageMerger} -i "${partial}" -j "${merged}" -o "${tmp}" 2> /dev/null > /dev/null || warning "error while calling ${txtImageMerger}"
+ echo "${header}" > "${merged}"
+ grep -v '## Merge' "${tmp}" >> "${merged}"
+ rm "${tmp}"
+ done
+ end_bar
+ echo " ${indent}merged ${count} files"
}
hdrImageMerger="clitkImageArithm"
test -x "./clitkImageArithm" && hdrImageMerger="./clitkImageArithm"
function merge_hdr_image {
-local merged="$1"
-local merged_bin="${merged%.*}.img"
-shift
-echo " ${indent}entering hdr image merger"
-echo " ${indent}merger is ${hdrImageMerger}"
-echo " ${indent}creating ${merged}"
-local count=0
-start_bar $#
-while test $# -gt 0
-do
- local partial="$1"
- local partial_bin="${partial%.*}.img"
+ local merged="$1"
+ local merged_bin="${merged%.*}.img"
shift
- let count++
-
- if test ! -f "${merged}"
- then
- update_bar ${count} "copying first partial result ${partial}"
- cp "${partial}" "${merged}"
- cp "${partial_bin}" "${merged_bin}"
- continue
- fi
+ echo " ${indent}entering hdr image merger"
+ echo " ${indent}merger is ${hdrImageMerger}"
+ echo " ${indent}creating ${merged}"
+ local count=0
+ start_bar $#
+ while test $# -gt 0
+ do
+ local partial="$1"
+ local partial_bin="${partial%.*}.img"
+ shift
+ let count++
+
+ if test ! -f "${merged}"
+ then
+ update_bar ${count} "copying first partial result ${partial}"
+ cp "${partial}" "${merged}"
+ cp "${partial_bin}" "${merged_bin}"
+ continue
+ fi
- update_bar ${count} "adding ${partial}"
- ${hdrImageMerger} -t 0 -i "${partial}" -j "${merged}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${hdrImageMerger}"
-done
-end_bar
-echo " ${indent}merged ${count} files"
+ update_bar ${count} "adding ${partial}"
+ ${hdrImageMerger} -t 0 -i "${partial}" -j "${merged}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${hdrImageMerger}"
+ done
+ end_bar
+ echo " ${indent}merged ${count} files"
}
mhdImageMerger="clitkImageArithm"
test -x "./clitkImageArithm" && mhdImageMerger="./clitkImageArithm"
function merge_mhd_image {
-local merged="$1"
-local merged_bin="${merged%.*}.raw"
-shift
-echo " ${indent}entering mhd image merger"
-echo " ${indent}merger is ${mhdImageMerger}"
-echo " ${indent}creating ${merged}"
-local count=0
-start_bar $#
-while test $# -gt 0
-do
- local partial="$1"
- local partial_bin="$(dirname "${partial}")/$(awk -F' = ' '/ElementDataFile/ { print $2; }' "${partial}")"
+ local merged="$1"
+ local merged_bin="${merged%.*}.raw"
shift
- let count++
-
- if test ! -f "${merged}"
- then
- update_bar ${count} "copying first partial result ${partial}"
- cp "${partial}" "${merged}"
- cp "${partial_bin}" "${merged_bin%.*}.${partial_bin##*.}"
- continue
- fi
+ echo " ${indent}entering mhd image merger"
+ echo " ${indent}merger is ${mhdImageMerger}"
+ echo " ${indent}creating ${merged}"
+ local count=0
+ start_bar $#
+ while test $# -gt 0
+ do
+ local partial="$1"
+ local partial_bin="$(dirname "${partial}")/$(awk -F' = ' '/ElementDataFile/ { print $2; }' "${partial}")"
+ shift
+ let count++
+
+ if test ! -f "${merged}"
+ then
+ update_bar ${count} "copying first partial result ${partial}"
+ cp "${partial}" "${merged}"
+ cp "${partial_bin}" "${merged_bin%.*}.${partial_bin##*.}"
+ continue
+ fi
- update_bar ${count} "adding ${partial}"
- ${mhdImageMerger} -t 0 -i "${partial}" -j "${merged}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${mhdImageMerger}"
- mv "${merged_bin}" "${merged_bin%.*}.${partial_bin##*.}"
- sed -i "s/$(basename "${merged_bin}")/$(basename "${merged_bin%.*}.${partial_bin##*.}")/" "${merged}"
-done
-end_bar
-echo " ${indent}merged ${count} files"
+ update_bar ${count} "adding ${partial}"
+ ${mhdImageMerger} -t 0 -i "${partial}" -j "${merged}" -o "${merged}" 2> /dev/null > /dev/null || warning "error while calling ${mhdImageMerger}"
+ mv "${merged_bin}" "${merged_bin%.*}.${partial_bin##*.}"
+ sed -i "s/$(basename "${merged_bin}")/$(basename "${merged_bin%.*}.${partial_bin##*.}")/" "${merged}"
+ done
+ end_bar
+ echo " ${indent}merged ${count} files"
}
function merge_dispatcher {
error "unknown file type"
}
+function merge_dispatcher_uncertainty {
+ local indent=" ** "
+ local outputfile="${1:?"provide output filename"}"
+
+ local partialoutputfiles="$(find "${rundir}" -mindepth 2 -type f -name "${outputfile}")"
+ local nboutputfiles="$(echo "${partialoutputfiles}" | wc -l)"
+ if test ${nboutputdirs} -ne ${nboutputfiles}
+ then
+ warning "missing files"
+ return
+ fi
+
+ local firstpartialoutputfile="$(echo "${partialoutputfiles}" | head -n 1)"
+ local firstpartialoutputextension="${firstpartialoutputfile##*.}"
+
+ if [[ "${firstpartialoutputfile}" == *Uncertainty* ]]
+ then
+ if test "${firstpartialoutputextension}" == "mhd"
+ then
+ echo "${indent}Uncertainty file found: ${firstpartialoutputfile}"
+ ## search for sum
+ local mergedfile="${outputdir}/$(basename "${firstpartialoutputfile}")"
+ summed_merged_file=${mergedfile//-Uncertainty/}
+ if [ ! -f ${summed_merged_file} ];
+ then
+ warning "${summed_merged_file} does not exist. Error, no uncertainty computed"
+ return;
+ fi
+ echo "${indent}${summed_merged_file} found"
+ ## search for Squared
+ squared_merged_file=${mergedfile//-Uncertainty/-Squared}
+ if [ ! -f ${squared_merged_file} ];
+ then
+ warning "${squared_merged_file} does not exist. Error, no uncertainty computed"
+ return;
+ fi
+ echo "${indent}${squared_merged_file} found"
+ ## search for NumberOfEvent
+ totalEvents=0;
+ for outputfile in $(find "${rundir}" -regextype 'posix-extended' -type f -regex "${rundir}/output.*\.(hdr|mhd|root|txt)" | awk -F '/' '{ print $NF; }' | sort | uniq)
+ do
+ #echo $outputfile
+ if grep -q 'NumberOfEvent' "${outputdir}/${outputfile}"
+ then
+ totalEvents="$(grep "NumberOfEvents" "${outputdir}/${outputfile}" | cut -d' ' -f4)"
+ echo "${indent}Find the NumberOfEvent in $outputfile: ${totalEvents}"
+ fi
+ done
+
+ if test ${totalEvents} -gt 0
+ then
+ uncerImageMerger="clitkImageUncertainty"
+ test -x "./clitkImageUncertainty" && uncerImageMerger="./clitkImageUncertainty"
+ ${uncerImageMerger} -i ${summed_merged_file} -s ${squared_merged_file} -o ${mergedfile} -n ${totalEvents}
+ else
+ warning "${totalEvents} not positive. A at least one stat file (SimulationStatisticActor) must be provided. Error, no uncertainty computed"
+ return;
+ fi
+ fi
+ fi
+
+}
+
echo "!!!! this is $0 v0.3k !!!!"
rundir="${1?"provide run dir"}"
fi
outputdir="$(basename "${outputdir}")"
echo "output dir is ${outputdir}"
+
test -d "${outputdir}" && rm -r "${outputdir}"
mkdir "${outputdir}"
merge_dispatcher "${outputfile}"
done
+echo ""
+echo "Merging done. Special case for statistical uncertainty"
+for outputfile in $(find "${outputdir}" -regextype 'posix-extended' -type f -regex "${outputdir}/output.*\.(hdr|mhd|root|txt)" | awk -F '/' '{ print $NF; }' | sort | uniq)
+do
+ merge_dispatcher_uncertainty "${outputfile}"
+done
+
if [ -f "${rundir}/params.txt" ]
then
- echo "copying params file"
- cp "${rundir}/params.txt" "${outputdir}/params.txt"
+ echo "copying params file"
+ cp "${rundir}/params.txt" "${outputdir}/params.txt"
fi
echo "these was ${warning_count} warning(s)"
clitkExceptionObject.cxx
clitkFilterBase.cxx
clitkMemoryUsage.cxx
+ clitkMatrix.cxx
vvImage.cxx
vvImageReader.cxx
vvImageWriter.cxx
# include <stdint.h>
#endif
+#define VTK_EXCLUDE_STRSTREAM_HEADERS
+
//--------------------------------------------------------------------
namespace clitk {
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+
+#ifndef clitkElastix_h
+#define clitkElastix_h
+
+#include <itkEuler3DTransform.h>
+#include <iterator>
+
+//--------------------------------------------------------------------
+namespace clitk {
+
+//-------------------------------------------------------------------
+bool
+GetElastixValueFromTag(std::ifstream & is,
+ std::string tag,
+ std::string & value)
+{
+ std::string line;
+ is.seekg (0, is.beg);
+ while(std::getline(is, line)) {
+ unsigned pos = line.find(tag);
+ if (pos<line.size()) {
+ value=line.substr(pos+tag.size(),line.size()-2);// remove last ')'
+ value.erase (std::remove (value.begin(), value.end(), '"'), value.end());
+ value.erase (std::remove (value.begin(), value.end(), ')'), value.end());
+ return true;
+ }
+ }
+ return false;
+}
+//-------------------------------------------------------------------
+
+
+//-------------------------------------------------------------------
+void
+GetValuesFromValue(const std::string & s,
+ std::vector<std::string> & values)
+{
+ std::stringstream strstr(s);
+ std::istream_iterator<std::string> it(strstr);
+ std::istream_iterator<std::string> end;
+ std::vector<std::string> results(it, end);
+ values.clear();
+ values.resize(results.size());
+ for(uint i=0; i<results.size(); i++) values[i] = results[i];
+}
+//-------------------------------------------------------------------
+
+
+//-------------------------------------------------------------------
+template<unsigned int Dimension>
+typename itk::Matrix<double, Dimension+1, Dimension+1>
+createMatrixFromElastixFile(std::vector<std::string> & filename, bool verbose=true) {
+ if (Dimension != 3) {
+ FATAL("Only 3D yet" << std::endl);
+ }
+ typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
+
+ itk::Euler3DTransform<double>::Pointer mat = itk::Euler3DTransform<double>::New();
+ itk::Euler3DTransform<double>::Pointer previous;
+ for(uint j=0; j<filename.size(); j++) {
+
+ // Open file
+ if (verbose) std::cout << "Read elastix parameters in " << filename[j] << std::endl;
+ std::ifstream is;
+ clitk::openFileForReading(is, filename[j]);
+
+ // Check Transform
+ std::string s;
+ bool b = GetElastixValueFromTag(is, "Transform ", s);
+ if (!b) {
+ FATAL("Error must read 'Transform' in " << filename[j] << std::endl);
+ }
+ if (s != "EulerTransform") {
+ FATAL("Sorry only 'EulerTransform'" << std::endl);
+ }
+
+ // FIXME check
+ // (InitialTransformParametersFilename[j] "NoInitialTransform")
+
+ // Get CenterOfRotationPoint
+ GetElastixValueFromTag(is, "CenterOfRotationPoint ", s); // space is needed
+ if (!b) {
+ FATAL("Error must read 'CenterOfRotationPoint' in " << filename[j] << std::endl);
+ }
+ std::vector<std::string> cor;
+ GetValuesFromValue(s, cor);
+ itk::Euler3DTransform<double>::CenterType c;
+ for(uint i=0; i<3; i++)
+ c[i] = atof(cor[i].c_str());
+ mat->SetCenter(c);
+
+ // Get Transformparameters
+ GetElastixValueFromTag(is, "ComputeZYX ", s); // space is needed
+ mat->SetComputeZYX( s==std::string("true") );
+
+ // Get Transformparameters
+ GetElastixValueFromTag(is, "TransformParameters ", s); // space is needed
+ if (!b) {
+ FATAL("Error must read 'TransformParameters' in " << filename[j] << std::endl);
+ }
+ std::vector<std::string> results;
+ GetValuesFromValue(s, results);
+
+ // construct a stream from the string
+ itk::Euler3DTransform<double>::ParametersType p;
+ p.SetSize(6);
+ for(uint i=0; i<3; i++)
+ p[i] = atof(results[i].c_str()); // Rotation
+ for(uint i=0; i<3; i++)
+ p[i+3] = atof(results[i+3].c_str()); // Translation
+ mat->SetParameters(p);
+
+ if (verbose) {
+ std::cout << "Rotation (deg) : " << rad2deg(p[0]) << " " << rad2deg(p[1]) << " " << rad2deg(p[2]) << std::endl;
+ std::cout << "Center of rot (phy) : " << c[0] << " " << c[1] << " " << c[2] << std::endl;
+ std::cout << "Translation (phy) : " << p[3] << " " << p[4] << " " << p[5] << std::endl;
+ }
+
+ // Compose with previous if needed
+ if (j!=0) {
+ mat->Compose(previous);
+ if (verbose) {
+ std::cout << "Composed rotation (deg) : " << rad2deg(mat->GetAngleX()) << " " << rad2deg(mat->GetAngleY()) << " " << rad2deg(mat->GetAngleZ()) << std::endl;
+ std::cout << "Composed center of rot (phy) : " << mat->GetCenter() << std::endl;
+ std::cout << "Compsoed translation (phy) : " << mat->GetTranslation() << std::endl;
+ }
+ }
+ // previous = mat->Clone(); // ITK4
+ previous = itk::Euler3DTransform<double>::New();
+ previous->SetParameters(mat->GetParameters());
+ previous->SetCenter(c);
+ previous->SetComputeZYX(mat->GetComputeZYX());
+ }
+
+ mat = previous;
+ for(uint i=0; i<3; i++)
+ for(uint j=0; j<3; j++)
+ matrix[i][j] = mat->GetMatrix()[i][j];
+ // Offset is -Rc + t + c
+ matrix[0][3] = mat->GetOffset()[0];
+ matrix[1][3] = mat->GetOffset()[1];
+ matrix[2][3] = mat->GetOffset()[2];
+ matrix[3][3] = 1;
+
+ return matrix;
+}
+}
+//-------------------------------------------------------------------
+
+#endif
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+
+#include "clitkMatrix.h"
+
+//--------------------------------------------------------------------
+namespace clitk {
+
+//-------------------------------------------------------------------
+std::string
+Get4x4MatrixDoubleAsString(vtkMatrix4x4 *matrix,
+ const int precision)
+{
+ std::ostringstream strmatrix;
+
+ // Figure out the number of digits of the integer part of the largest absolute value
+ // for each column
+ unsigned width[4];
+ for (unsigned int j = 0; j < 4; j++){
+ double absmax = 0.;
+ for (unsigned int i = 0; i < 4; i++)
+ absmax = std::max(absmax, vnl_math_abs(matrix->GetElement(i, j)));
+ unsigned ndigits = (unsigned)std::max(0.,std::log10(absmax))+1;
+ width[j] = precision+ndigits+3;
+ }
+
+ // Output with correct width, aligned to the right
+ for (unsigned int i = 0; i < 4; i++) {
+ for (unsigned int j = 0; j < 4; j++) {
+ strmatrix.setf(ios::fixed,ios::floatfield);
+ strmatrix.precision(precision);
+ strmatrix.fill(' ');
+ strmatrix.width(width[j]);
+ strmatrix << std::right << matrix->GetElement(i, j);
+ }
+ strmatrix << std::endl;
+ }
+ std::string result = strmatrix.str().c_str();
+ return result;
+}
+//-------------------------------------------------------------------
+
+//-------------------------------------------------------------------
+std::string
+Get4x4MatrixDoubleAsString(itk::Matrix<double, 4, 4> m,
+ const int precision)
+{
+ vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
+ for (unsigned int j = 0; j < 4; j++)
+ for (unsigned int i = 0; i < 4; i++)
+ matrix->SetElement(j,i,m[j][i]);
+ return Get4x4MatrixDoubleAsString(matrix, precision);
+}
+//-------------------------------------------------------------------
+}
+//-------------------------------------------------------------------
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+
+#ifndef clitkMatrix_h
+#define clitkMatrix_h
+
+#include <itkMatrix.h>
+#define VTK_EXCLUDE_STRSTREAM_HEADERS
+#include <vtkMatrix4x4.h>
+#include <vtkSmartPointer.h>
+
+//--------------------------------------------------------------------
+namespace clitk {
+std::string
+Get4x4MatrixDoubleAsString(vtkMatrix4x4 *matrix,
+ const int precision=3);
+
+std::string
+Get4x4MatrixDoubleAsString(itk::Matrix<double, 4, 4> m,
+ const int precision=3);
+}
+//-------------------------------------------------------------------
+
+#endif
#include "itkPoint.h"
#include "clitkImageCommon.h"
#include "clitkCommon.h"
+#define VTK_EXCLUDE_STRSTREAM_HEADERS
#include <vtkMatrix4x4.h>
#include <vtkSmartPointer.h>
#include <itkObjectFactory.h>
#include <itkProcessObject.h>
+#define VTK_EXCLUDE_STRSTREAM_HEADERS
#include <vtkSmartPointer.h>
#include <vtkTransform.h>
#include "vvImageReader.h"
#include "vvImageReader.txx"
#include "clitkTransformUtilities.h"
+#include "clitkElastix.h"
//------------------------------------------------------------------------------
vvImageReader::vvImageReader()
{
std::string filename(mInputFilenames[0]);
std::string ext(itksys::SystemTools::GetFilenameLastExtension(filename));
+
+ // Try a ".mat" extension
if (ext.length() > 0) {
size_t pos = filename.rfind(ext);
filename.replace(pos, ext.length(), ".mat");
}
else
filename += ".mat";
-
std::ifstream f(filename.c_str());
+ itk::Matrix<double, 4, 4> itkMat;
+ bool itkMatRead = false;
if(f.is_open()) {
- f.close();
+ itkMatRead = true;
- itk::Matrix<double, 4, 4> itkMat;
itkMat.SetIdentity();
try {
itkMat = clitk::ReadMatrix3D(filename);
catch (itk::ExceptionObject & err) {
itkWarningMacro(<< "Found " << filename
<< " but this is not a 4x4 matrix so it is ignored.");
+ itkMatRead = false;
}
+ }
+ f.close();
+
+ // Try a ".elx" extension
+ filename = mInputFilenames[0];
+ if (ext.length() > 0) {
+ size_t pos = filename.rfind(ext);
+ filename.replace(pos, ext.length(), ".elx");
+ }
+ else
+ filename += ".elx";
+ f.open(filename.c_str());
+ if(!itkMatRead && f.is_open()) {
+ itkMatRead = true;
+ std::vector<std::string> l;
+ l.push_back(filename);
+ itkMat = clitk::createMatrixFromElastixFile<3>(l, true);
+ }
+ f.close();
+ if(itkMatRead) {
vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
matrix->Identity();
for(int j=0; j<4; j++)
//for image sequences, apply the transform to each images of the sequence
if (mImage->IsTimeSequence())
{
- for (unsigned i = 1 ; i<mImage->GetTransform().size() ; i++)
- {
- mImage->GetTransform()[i]->PreMultiply();
- mImage->GetTransform()[i]->Concatenate(matrix);
- mImage->GetTransform()[i]->Update();
- }
+ for (unsigned i = 1 ; i<mImage->GetTransform().size() ; i++)
+ {
+ mImage->GetTransform()[i]->PreMultiply();
+ mImage->GetTransform()[i]->Concatenate(matrix);
+ mImage->GetTransform()[i]->Update();
+ }
}
}
/*=========================================================================
Program: vv http://www.creatis.insa-lyon.fr/rio/vv
- Authors belong to:
+ Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
#include "clitkFilterBase.h"
#include "clitkCropLikeImageFilter.h"
#include "clitkSegmentationUtils.h"
-// #include "itkLabelOverlapMeasuresImageFilter.h"
// itk
#include <itkImageToImageFilter.h>
#include <itkLabelStatisticsImageFilter.h>
+#include <itkLabelOverlapMeasuresImageFilter.h>
#include <iomanip>
namespace clitk {
-
+
//--------------------------------------------------------------------
/*
Analyze the relative position of a Target (mask image) according
to some Object, in a given Support. Indicate the main important
- position of this Target according the Object.
+ position of this Target according the Object.
*/
//--------------------------------------------------------------------
-
+
template <class ImageType>
class LabelImageOverlapMeasureFilter:
public virtual FilterBase,
typedef LabelImageOverlapMeasureFilter<ImageType> Self;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
-
+
/** Method for creation through the object factory. */
itkNewMacro(Self);
-
+
/** Run-time type information (and related methods). */
itkTypeMacro(LabelImageOverlapMeasureFilter, ImageToImageFilter);
/** Some convenient typedefs. */
typedef typename ImageType::ConstPointer ImageConstPointer;
typedef typename ImageType::Pointer ImagePointer;
- typedef typename ImageType::RegionType RegionType;
+ typedef typename ImageType::RegionType RegionType;
typedef typename ImageType::PixelType PixelType;
typedef typename ImageType::SpacingType SpacingType;
typedef typename ImageType::SizeType SizeType;
typedef typename ImageType::IndexType IndexType;
typedef typename ImageType::PointType PointType;
-
+
/** ImageDimension constants */
itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension);
FILTERBASE_INIT;
-
+
// Options
itkGetConstMacro(BackgroundValue, PixelType);
- itkSetMacro(BackgroundValue, PixelType);
+ //itkSetMacro(BackgroundValue, PixelType);
itkGetConstMacro(Label1, PixelType);
itkSetMacro(Label1, PixelType);
itkGetConstMacro(Label2, PixelType);
itkSetMacro(Label2, PixelType);
+ itkSetMacro(VerboseFlag, bool);
+ itkGetConstMacro(VerboseFlag, bool);
+ itkBooleanMacro(VerboseFlag);
+
+ itkSetMacro(LongFlag, bool);
+ itkGetConstMacro(LongFlag, bool);
+ itkBooleanMacro(LongFlag);
+
// I dont want to verify inputs information
virtual void VerifyInputInformation() { }
-
+
protected:
LabelImageOverlapMeasureFilter();
virtual ~LabelImageOverlapMeasureFilter() {}
-
+
PixelType m_BackgroundValue;
PixelType m_Label1;
PixelType m_Label2;
ImagePointer m_Input1;
ImagePointer m_Input2;
+ bool m_VerboseFlag;
+ bool m_LongFlag;
virtual void GenerateOutputInformation();
virtual void GenerateInputRequestedRegion();
private:
LabelImageOverlapMeasureFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
-
+
}; // end class
//--------------------------------------------------------------------
/*=========================================================================
Program: vv http://www.creatis.insa-lyon.fr/rio/vv
- Authors belong to:
+ Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
this->SetNumberOfRequiredInputs(2);
SetLabel1(1);
SetLabel2(1);
- SetBackgroundValue(0);
+ m_BackgroundValue = 0;
+ SetVerboseFlag(false);
+ SetLongFlag(false);
}
//--------------------------------------------------------------------
//--------------------------------------------------------------------
template <class ImageType>
-void
+void
clitk::LabelImageOverlapMeasureFilter<ImageType>::
-GenerateOutputInformation()
-{
+GenerateOutputInformation()
+{
// DD("GenerateOutputInformation");
//ImagePointer input = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
// ImagePointer outputImage = this->GetOutput(0);
//--------------------------------------------------------------------
template <class ImageType>
-void
+void
clitk::LabelImageOverlapMeasureFilter<ImageType>::
-GenerateInputRequestedRegion()
+GenerateInputRequestedRegion()
{
// Call default
itk::ImageToImageFilter<ImageType, ImageType>::GenerateInputRequestedRegion();
//--------------------------------------------------------------------
template <class ImageType>
-void
+void
clitk::LabelImageOverlapMeasureFilter<ImageType>::
-GenerateData()
+GenerateData()
{
- // DD("GenerateData");
-
// Get input pointer
m_Input1 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(0));
m_Input2 = dynamic_cast<ImageType*>(itk::ProcessObject::GetInput(1));
// Resize like the union
ImagePointer input1 = clitk::ResizeImageLike<ImageType>(m_Input1, bbo, GetBackgroundValue());
ImagePointer input2 = clitk::ResizeImageLike<ImageType>(m_Input2, bbo, GetBackgroundValue());
- //DD(input1->GetLargestPossibleRegion());
- //DD(input2->GetLargestPossibleRegion());
-
- // Compute overlap image
- ImagePointer image_union = clitk::Clone<ImageType>(input1);
- ImagePointer image_intersection = clitk::Clone<ImageType>(input1);
- clitk::Or<ImageType>(image_union, input2, GetBackgroundValue());
- clitk::And<ImageType>(image_intersection, input2, GetBackgroundValue());
-
- ImagePointer image_1NotIn2 = clitk::Clone<ImageType>(input1);
- clitk::AndNot<ImageType>(image_1NotIn2, input2, GetBackgroundValue());
-
- ImagePointer image_2NotIn1 = clitk::Clone<ImageType>(input2);
- clitk::AndNot<ImageType>(image_2NotIn1, input1, GetBackgroundValue());
-
- //writeImage<ImageType>(image_union, "union.mha");
- //writeImage<ImageType>(image_intersection, "intersection.mha");
- //writeImage<ImageType>(image_1NotIn2, "image_1NotIn2.mha");
- //writeImage<ImageType>(image_2NotIn1, "image_2NotIn1.mha");
-
- // Compute size
- typedef itk::LabelStatisticsImageFilter<ImageType, ImageType> StatFilterType;
- typename StatFilterType::Pointer statFilter = StatFilterType::New();
- statFilter->SetInput(image_union);
- statFilter->SetLabelInput(image_union);
- statFilter->Update();
- int u = statFilter->GetCount(GetLabel1());
-
- statFilter->SetInput(image_intersection);
- statFilter->SetLabelInput(image_intersection);
- statFilter->Update();
- int inter = statFilter->GetCount(GetLabel1());
-
- statFilter->SetInput(m_Input1);
- statFilter->SetLabelInput(m_Input1);
- statFilter->Update();
- int in1 = statFilter->GetCount(GetLabel1());
-
- statFilter->SetInput(m_Input2);
- statFilter->SetLabelInput(m_Input2);
- statFilter->Update();
- int in2 = statFilter->GetCount(GetLabel1());
-
- statFilter->SetInput(image_1NotIn2);
- statFilter->SetLabelInput(image_1NotIn2);
- statFilter->Update();
- int l1notIn2 = statFilter->GetCount(GetLabel1());
-
- statFilter->SetInput(image_2NotIn1);
- statFilter->SetLabelInput(image_2NotIn1);
- statFilter->Update();
- int l2notIn1 = statFilter->GetCount(GetLabel1());
-
- double dice = 2.0*(double)inter/(double)(in1+in2);
+ /*
+ if (GetVerboseFlag()) {
+ std::cout << "Resize images to the union of bounding boxes: "
+ << input1->GetLargestPossibleRegion().GetSize() << std::endl;
+ }
+ */
+
+ /*
int width = 6;
std::cout << std::setw(width) << in1 << " "
<< std::setw(width) << in2 << " "
<< std::setw(width) << l1notIn2 << " "
<< std::setw(width) << l2notIn1 << " "
<< std::setw(width) << dice << " "; //std::endl;
+ */
+
+ // Compute overlap, dice
+ typedef itk::LabelOverlapMeasuresImageFilter<ImageType> FilterType;
+ typename FilterType::Pointer diceFilter = FilterType::New();
+ diceFilter->SetSourceImage(input1);
+ diceFilter->SetTargetImage(input2);
+ diceFilter->Update();
+
+ // Display information
+ if (!GetLongFlag()) {
+ if (GetVerboseFlag()) {
+ std::cout << "Dice : ";
+ }
+ std::cout << diceFilter->GetDiceCoefficient() << std::endl;
+ }
+ else { // long options
+ if (GetVerboseFlag()) {
+ std::cout << "Dice Jaccard Source Target Inter Union SrComp TarComp" << std::endl;
+ }
+ typename FilterType::MapType m = diceFilter->GetLabelSetMeasures();
+ int width = 6;
+ std::cout << std::setw(width) << diceFilter->GetDiceCoefficient() << " "
+ << std::setw(width) << diceFilter->GetJaccardCoefficient() << " "
+ << std::setw(width) << m[m_Label1].m_Source << " "
+ << std::setw(width) << m[m_Label2].m_Target << " "
+ << std::setw(width) << m[m_Label1].m_Intersection << " "
+ << std::setw(width) << m[m_Label1].m_Union << " "
+ << std::setw(width) << m[m_Label1].m_SourceComplement << " "
+ << std::setw(width) << m[m_Label2].m_TargetComplement << " "
+ << std::endl;
+ }
+
}
//--------------------------------------------------------------------
-
std::cout << "LastDimIsTime = " << m_LastDimensionIsTime << std::endl;
}
+ // Compute origin based on image corner
+ typename FilterType::OriginPointType origin = input->GetOrigin();
+ for(unsigned int i=0; i<OutputImageType::ImageDimension; i++) {
+ origin[i] -= 0.5 * input->GetSpacing()[i];
+ origin[i] += 0.5 * m_OutputSpacing[i];
+ }
+
// Instance of the transform object to be passed to the resample
// filter. By default, identity transform is applied
filter->SetTransform(m_Transform);
filter->SetSize(m_OutputSize);
filter->SetOutputSpacing(m_OutputSpacing);
- filter->SetOutputOrigin(input->GetOrigin());
+ filter->SetOutputOrigin(origin);
filter->SetDefaultPixelValue(m_DefaultPixelValue);
filter->SetNumberOfThreads(this->GetNumberOfThreads());
filter->SetOutputDirection(input->GetDirection()); // <-- NEEDED if we want to keep orientation (in case of PermutAxes for example)
+++ /dev/null
-/*=========================================================================
-
- Program: Insight Segmentation & Registration Toolkit
- Module: $RCSfile: itkLabelOverlapMeasuresImageFilter.h,v $
- Language: C++
- Date: $Date: $
- Version: $Revision: $
-
- Copyright (c) Insight Software Consortium. All rights reserved.
- See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-#ifndef __itkLabelOverlapMeasuresImageFilter_h
-#define __itkLabelOverlapMeasuresImageFilter_h
-
-#include "itkImageToImageFilter.h"
-#include "itkFastMutexLock.h"
-#include "itkNumericTraits.h"
-
-//#include "itk_hash_map.h"
-#include "itksys/hash_map.hxx"
-
-namespace itk {
-
-/** \class LabelOverlapMeasuresImageFilter
- * \brief Computes overlap measures between the set same set of labels of
- * pixels of two images. Background is assumed to be 0.
- *
- * \sa LabelOverlapMeasuresImageFilter
- *
- * \ingroup MultiThreaded
- */
-template<class TLabelImage>
-class ITK_EXPORT LabelOverlapMeasuresImageFilter :
- public ImageToImageFilter<TLabelImage, TLabelImage>
-{
-public:
- /** Standard Self typedef */
- typedef LabelOverlapMeasuresImageFilter Self;
- typedef ImageToImageFilter<TLabelImage,TLabelImage> Superclass;
- typedef SmartPointer<Self> Pointer;
- typedef SmartPointer<const Self> ConstPointer;
-
- /** Method for creation through the object factory. */
- itkNewMacro( Self );
-
- /** Runtime information support. */
- itkTypeMacro( LabelOverlapMeasuresImageFilter, ImageToImageFilter );
-
- virtual void VerifyInputInformation() { }
-
-
- /** Image related typedefs. */
- typedef TLabelImage LabelImageType;
- typedef typename TLabelImage::Pointer LabelImagePointer;
- typedef typename TLabelImage::ConstPointer LabelImageConstPointer;
-
- typedef typename TLabelImage::RegionType RegionType;
- typedef typename TLabelImage::SizeType SizeType;
- typedef typename TLabelImage::IndexType IndexType;
-
- typedef typename TLabelImage::PixelType LabelType;
-
- /** Type to use form computations. */
- typedef typename NumericTraits<LabelType>::RealType RealType;
-
- /** \class LabelLabelOverlapMeasuress
- * \brief Metrics stored per label */
- class LabelSetMeasures
- {
- public:
- // default constructor
- LabelSetMeasures()
- {
- m_Source = 0;
- m_Target = 0;
- m_Union = 0;
- m_Intersection = 0;
- m_SourceComplement = 0;
- m_TargetComplement = 0;
- }
-
- // added for completeness
- LabelSetMeasures& operator=( const LabelSetMeasures& l )
- {
- m_Source = l.m_Source;
- m_Target = l.m_Target;
- m_Union = l.m_Union;
- m_Intersection = l.m_Intersection;
- m_SourceComplement = l.m_SourceComplement;
- m_TargetComplement = l.m_TargetComplement;
- }
-
- unsigned long m_Source;
- unsigned long m_Target;
- unsigned long m_Union;
- unsigned long m_Intersection;
- unsigned long m_SourceComplement;
- unsigned long m_TargetComplement;
- };
-
- /** Type of the map used to store data per label */
- typedef itksys::hash_map<LabelType, LabelSetMeasures> MapType;
- typedef typename MapType::iterator MapIterator;
- typedef typename MapType::const_iterator MapConstIterator;
-
- /** Image related typedefs. */
- itkStaticConstMacro( ImageDimension, unsigned int,
- TLabelImage::ImageDimension );
-
- /** Set the source image. */
- void SetSourceImage( const LabelImageType * image )
- { this->SetNthInput( 0, const_cast<LabelImageType *>( image ) ); }
-
- /** Set the target image. */
- void SetTargetImage( const LabelImageType * image )
- { this->SetNthInput( 1, const_cast<LabelImageType *>( image ) ); }
-
- /** Get the source image. */
- const LabelImageType * GetSourceImage( void )
- { return this->GetInput( 0 ); }
-
- /** Get the target image. */
- const LabelImageType * GetTargetImage( void )
- { return this->GetInput( 1 ); }
-
- /** Get the label set measures */
- MapType GetLabelSetMeasures()
- { return this->m_LabelSetMeasures; }
-
- /**
- * tric overlap measures
- */
- /** measures over all labels */
- RealType GetTotalOverlap();
- RealType GetUnionOverlap();
- RealType GetMeanOverlap();
- RealType GetVolumeSimilarity();
- RealType GetFalseNegativeError();
- RealType GetFalsePositiveError();
- /** measures over individual labels */
- RealType GetTargetOverlap( LabelType );
- RealType GetUnionOverlap( LabelType );
- RealType GetMeanOverlap( LabelType );
- RealType GetVolumeSimilarity( LabelType );
- RealType GetFalseNegativeError( LabelType );
- RealType GetFalsePositiveError( LabelType );
- /** alternative names */
- RealType GetJaccardCoefficient()
- { return this->GetUnionOverlap(); }
- RealType GetJaccardCoefficient( LabelType label )
- { return this->GetUnionOverlap( label ); }
- RealType GetDiceCoefficient()
- { return this->GetMeanOverlap(); }
- RealType GetDiceCoefficient( LabelType label )
- { return this->GetMeanOverlap( label ); }
-
-
-#ifdef ITK_USE_CONCEPT_CHECKING
- /** Begin concept checking */
- itkConceptMacro( Input1HasNumericTraitsCheck,
- ( Concept::HasNumericTraits<LabelType> ) );
- /** End concept checking */
-#endif
-
-protected:
- LabelOverlapMeasuresImageFilter();
- ~LabelOverlapMeasuresImageFilter(){};
- void PrintSelf( std::ostream& os, Indent indent ) const;
-
- /**
- * Pass the input through unmodified. Do this by setting the output to the
- * source this by setting the output to the source image in the
- * AllocateOutputs() method.
- */
- void AllocateOutputs();
-
- void BeforeThreadedGenerateData();
-
- void AfterThreadedGenerateData();
-
- /** Multi-thread version GenerateData. */
- void ThreadedGenerateData( const RegionType&, int );
-
- // Override since the filter needs all the data for the algorithm
- void GenerateInputRequestedRegion();
-
- // Override since the filter produces all of its output
- void EnlargeOutputRequestedRegion( DataObject *data );
-
-private:
- LabelOverlapMeasuresImageFilter( const Self& ); //purposely not implemented
- void operator=( const Self& ); //purposely not implemented
-
- std::vector<MapType> m_LabelSetMeasuresPerThread;
- MapType m_LabelSetMeasures;
-
- SimpleFastMutexLock m_Mutex;
-
-}; // end of class
-
-} // end namespace itk
-
-#ifndef ITK_MANUAL_INSTANTIATION
-#include "itkLabelOverlapMeasuresImageFilter.txx"
-#endif
-
-#endif
+++ /dev/null
-/*=========================================================================
-
- Program: Insight Segmentation & Registration Toolkit
- Module: $RCSfile: itkLabelOverlapMeasuresImageFilter.txx,v $
- Language: C++
- Date: $Date: $
- Version: $Revision: $
-
- Copyright (c) Insight Software Consortium. All rights reserved.
- See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-
- This software is distributed WITHOUT ANY WARRANTY; without even
- the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
- PURPOSE. See the above copyright notices for more information.
-
-=========================================================================*/
-#ifndef _itkLabelOverlapMeasuresImageFilter_txx
-#define _itkLabelOverlapMeasuresImageFilter_txx
-
-#include "itkLabelOverlapMeasuresImageFilter.h"
-
-#include "itkImageRegionConstIterator.h"
-#include "itkProgressReporter.h"
-
-namespace itk {
-
-#if defined(__GNUC__) && (__GNUC__ <= 2) //NOTE: This class needs a mutex for gnu 2.95
-/** Used for mutex locking */
-#define LOCK_HASHMAP this->m_Mutex.Lock()
-#define UNLOCK_HASHMAP this->m_Mutex.Unlock()
-#else
-#define LOCK_HASHMAP
-#define UNLOCK_HASHMAP
-#endif
-
-template<class TLabelImage>
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::LabelOverlapMeasuresImageFilter()
-{
- // this filter requires two input images
- this->SetNumberOfRequiredInputs( 2 );
-}
-
-template<class TLabelImage>
-void
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GenerateInputRequestedRegion()
-{
- Superclass::GenerateInputRequestedRegion();
- if( this->GetSourceImage() )
- {
- LabelImagePointer source = const_cast
- <LabelImageType *>( this->GetSourceImage() );
- source->SetRequestedRegionToLargestPossibleRegion();
- }
- if( this->GetTargetImage() )
- {
- LabelImagePointer target = const_cast
- <LabelImageType *>( this->GetTargetImage() );
- target->SetRequestedRegionToLargestPossibleRegion();
- }
-}
-
-template<class TLabelImage>
-void
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::EnlargeOutputRequestedRegion( DataObject *data )
-{
- Superclass::EnlargeOutputRequestedRegion( data );
- data->SetRequestedRegionToLargestPossibleRegion();
-}
-
-
-template<class TLabelImage>
-void
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::AllocateOutputs()
-{
- // Pass the source through as the output
- LabelImagePointer image =
- const_cast<TLabelImage *>( this->GetSourceImage() );
- this->SetNthOutput( 0, image );
-
- // Nothing that needs to be allocated for the remaining outputs
-}
-
-template<class TLabelImage>
-void
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::BeforeThreadedGenerateData()
-{
- int numberOfThreads = this->GetNumberOfThreads();
-
- // Resize the thread temporaries
- this->m_LabelSetMeasuresPerThread.resize( numberOfThreads );
-
- // Initialize the temporaries
- for( int n = 0; n < numberOfThreads; n++ )
- {
- this->m_LabelSetMeasuresPerThread[n].clear();
- }
-
- // Initialize the final map
- this->m_LabelSetMeasures.clear();
-}
-
-template<class TLabelImage>
-void
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::AfterThreadedGenerateData()
-{
- // Run through the map for each thread and accumulate the set measures.
- for( int n = 0; n < this->GetNumberOfThreads(); n++ )
- {
- // iterate over the map for this thread
- for( MapConstIterator threadIt = this->m_LabelSetMeasuresPerThread[n].begin();
- threadIt != this->m_LabelSetMeasuresPerThread[n].end();
- ++threadIt )
- {
- // does this label exist in the cumulative stucture yet?
- MapIterator mapIt = this->m_LabelSetMeasures.find( ( *threadIt ).first );
- if( mapIt == this->m_LabelSetMeasures.end() )
- {
- // create a new entry
- typedef typename MapType::value_type MapValueType;
- mapIt = this->m_LabelSetMeasures.insert( MapValueType(
- (*threadIt).first, LabelSetMeasures() ) ).first;
- }
-
- // accumulate the information from this thread
- (*mapIt).second.m_Source += (*threadIt).second.m_Source;
- (*mapIt).second.m_Target += (*threadIt).second.m_Target;
- (*mapIt).second.m_Union += (*threadIt).second.m_Union;
- (*mapIt).second.m_Intersection +=
- (*threadIt).second.m_Intersection;
- (*mapIt).second.m_SourceComplement +=
- (*threadIt).second.m_SourceComplement;
- (*mapIt).second.m_TargetComplement +=
- (*threadIt).second.m_TargetComplement;
- } // end of thread map iterator loop
- } // end of thread loop
-}
-
-template<class TLabelImage>
-void
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::ThreadedGenerateData( const RegionType& outputRegionForThread,
- int threadId )
-{
- ImageRegionConstIterator<LabelImageType> ItS( this->GetSourceImage(),
- outputRegionForThread );
- ImageRegionConstIterator<LabelImageType> ItT( this->GetTargetImage(),
- outputRegionForThread );
-
- // support progress methods/callbacks
- ProgressReporter progress( this, threadId,
- 2*outputRegionForThread.GetNumberOfPixels() );
-
- for( ItS.GoToBegin(), ItT.GoToBegin(); !ItS.IsAtEnd(); ++ItS, ++ItT )
- {
- LabelType sourceLabel = ItS.Get();
- LabelType targetLabel = ItT.Get();
-
- // is the label already in this thread?
- MapIterator mapItS =
- this->m_LabelSetMeasuresPerThread[threadId].find( sourceLabel );
- MapIterator mapItT =
- this->m_LabelSetMeasuresPerThread[threadId].find( targetLabel );
-
- if( mapItS == this->m_LabelSetMeasuresPerThread[threadId].end() )
- {
- // create a new label set measures object
- typedef typename MapType::value_type MapValueType;
- LOCK_HASHMAP;
- mapItS = this->m_LabelSetMeasuresPerThread[threadId].insert(
- MapValueType( sourceLabel, LabelSetMeasures() ) ).first;
- UNLOCK_HASHMAP;
- }
-
- if( mapItT == this->m_LabelSetMeasuresPerThread[threadId].end() )
- {
- // create a new label set measures object
- typedef typename MapType::value_type MapValueType;
- LOCK_HASHMAP;
- mapItT = this->m_LabelSetMeasuresPerThread[threadId].insert(
- MapValueType( targetLabel, LabelSetMeasures() ) ).first;
- UNLOCK_HASHMAP;
- }
-
- (*mapItS).second.m_Source++;
- (*mapItT).second.m_Target++;
-
- if( sourceLabel == targetLabel )
- {
- (*mapItS).second.m_Intersection++;
- (*mapItS).second.m_Union++;
- }
- else
- {
- (*mapItS).second.m_Union++;
- (*mapItT).second.m_Union++;
-
- (*mapItS).second.m_SourceComplement++;
- (*mapItT).second.m_TargetComplement++;
- }
-
- progress.CompletedPixel();
- }
-}
-
-/**
- * measures
- */
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetTotalOverlap()
-{
- RealType numerator = 0.0;
- RealType denominator = 0.0;
- for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
- mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
- {
- // Do not include the background in the final value.
- if( (*mapIt).first == NumericTraits<LabelType>::Zero )
- {
- continue;
- }
- numerator += static_cast<RealType>( (*mapIt).second.m_Intersection );
- denominator += static_cast<RealType>( (*mapIt).second.m_Target );
- }
- return ( numerator / denominator );
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetTargetOverlap( LabelType label )
-{
- MapIterator mapIt = this->m_LabelSetMeasures.find( label );
- if( mapIt == this->m_LabelSetMeasures.end() )
- {
- itkWarningMacro( "Label " << label << " not found." );
- return 0.0;
- }
- RealType value =
- static_cast<RealType>( (*mapIt).second.m_Intersection ) /
- static_cast<RealType>( (*mapIt).second.m_Target );
- return value;
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetUnionOverlap()
-{
- RealType numerator = 0.0;
- RealType denominator = 0.0;
- for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
- mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
- {
- // Do not include the background in the final value.
- if( (*mapIt).first == NumericTraits<LabelType>::Zero )
- {
- continue;
- }
- numerator += static_cast<RealType>( (*mapIt).second.m_Intersection );
- denominator += static_cast<RealType>( (*mapIt).second.m_Union );
- }
- return ( numerator / denominator );
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetUnionOverlap( LabelType label )
-{
- MapIterator mapIt = this->m_LabelSetMeasures.find( label );
- if( mapIt == this->m_LabelSetMeasures.end() )
- {
- itkWarningMacro( "Label " << label << " not found." );
- return 0.0;
- }
- RealType value =
- static_cast<RealType>( (*mapIt).second.m_Intersection ) /
- static_cast<RealType>( (*mapIt).second.m_Union );
- return value;
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetMeanOverlap()
-{
- RealType uo = this->GetUnionOverlap();
- return ( 2.0 * uo / ( 1.0 + uo ) );
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetMeanOverlap( LabelType label )
-{
- RealType uo = this->GetUnionOverlap( label );
- return ( 2.0 * uo / ( 1.0 + uo ) );
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetVolumeSimilarity()
-{
- RealType numerator = 0.0;
- RealType denominator = 0.0;
- for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
- mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
- {
- // Do not include the background in the final value.
- if( (*mapIt).first == NumericTraits<LabelType>::Zero )
- {
- continue;
- }
- numerator += ( ( static_cast<RealType>( (*mapIt).second.m_Source ) -
- static_cast<RealType>( (*mapIt).second.m_Target ) ) );
- denominator += ( ( static_cast<RealType>( (*mapIt).second.m_Source ) +
- static_cast<RealType>( (*mapIt).second.m_Target ) ) );
- }
- return ( 2.0 * numerator / denominator );
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetVolumeSimilarity( LabelType label )
-{
- MapIterator mapIt = this->m_LabelSetMeasures.find( label );
- if( mapIt == this->m_LabelSetMeasures.end() )
- {
- itkWarningMacro( "Label " << label << " not found." );
- return 0.0;
- }
- RealType value = 2.0 *
- ( static_cast<RealType>( (*mapIt).second.m_Source ) -
- static_cast<RealType>( (*mapIt).second.m_Target ) ) /
- ( static_cast<RealType>( (*mapIt).second.m_Source ) +
- static_cast<RealType>( (*mapIt).second.m_Target ) );
- return value;
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetFalseNegativeError()
-{
- RealType numerator = 0.0;
- RealType denominator = 0.0;
- for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
- mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
- {
- // Do not include the background in the final value.
- if( (*mapIt).first == NumericTraits<LabelType>::Zero )
- {
- continue;
- }
- numerator += static_cast<RealType>( (*mapIt).second.m_TargetComplement );
- denominator += static_cast<RealType>( (*mapIt).second.m_Target );
- }
- return ( numerator / denominator );
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetFalseNegativeError( LabelType label )
-{
- MapIterator mapIt = this->m_LabelSetMeasures.find( label );
- if( mapIt == this->m_LabelSetMeasures.end() )
- {
- itkWarningMacro( "Label " << label << " not found." );
- return 0.0;
- }
- RealType value =
- static_cast<RealType>( (*mapIt).second.m_TargetComplement ) /
- static_cast<RealType>( (*mapIt).second.m_Target );
- return value;
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetFalsePositiveError()
-{
- RealType numerator = 0.0;
- RealType denominator = 0.0;
- for( MapIterator mapIt = this->m_LabelSetMeasures.begin();
- mapIt != this->m_LabelSetMeasures.end(); ++mapIt )
- {
- // Do not include the background in the final value.
- if( (*mapIt).first == NumericTraits<LabelType>::Zero )
- {
- continue;
- }
- numerator += static_cast<RealType>( (*mapIt).second.m_SourceComplement );
- denominator += static_cast<RealType>( (*mapIt).second.m_Source );
- }
- return ( numerator / denominator );
-}
-
-template<class TLabelImage>
-typename LabelOverlapMeasuresImageFilter<TLabelImage>::RealType
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::GetFalsePositiveError( LabelType label )
-{
- MapIterator mapIt = this->m_LabelSetMeasures.find( label );
- if( mapIt == this->m_LabelSetMeasures.end() )
- {
- itkWarningMacro( "Label " << label << " not found." );
- return 0.0;
- }
- RealType value =
- static_cast<RealType>( (*mapIt).second.m_SourceComplement ) /
- static_cast<RealType>( (*mapIt).second.m_Source );
- return value;
-}
-
-template<class TLabelImage>
-void
-LabelOverlapMeasuresImageFilter<TLabelImage>
-::PrintSelf( std::ostream& os, Indent indent ) const
-{
- Superclass::PrintSelf( os, indent );
-
-}
-
-
-}// end namespace itk
-#endif
SET(VTK_DIR ${build_prefix}/VTK)
#=========================================================
+#=========================================================
+# GDCM
+ ExternalProject_Add(
+ GDCM
+ SOURCE_DIR ${source_prefix}/gdcm
+ GIT_REPOSITORY git://git.code.sf.net/p/gdcm/gdcm
+ GIT_TAG v2.2.4
+ INSTALL_COMMAND ""
+ CMAKE_ARGS
+ -DCMAKE_C_FLAGS_DEBUG:STRING=${CMAKE_C_FLAGS_DEBUG}
+ -DCMAKE_C_FLAGS_MINSIZEREL:STRING=${CMAKE_C_FLAGS_MINSIZEREL}
+ -DCMAKE_C_FLAGS_RELEASE:STRING=${CMAKE_C_FLAGS_RELEASE}
+ -DCMAKE_C_FLAGS_RELWITHDEBINFO:STRING=${CMAKE_C_FLAGS_RELWITHDEBINFO}
+ -DCMAKE_CXX_FLAGS_DEBUG:STRING=${CMAKE_CXX_FLAGS_DEBUG}
+ -DCMAKE_CXX_FLAGS_MINSIZEREL:STRING=${CMAKE_CXX_FLAGS_MINSIZEREL}
+ -DCMAKE_CXX_FLAGS_RELEASE:STRING=${CMAKE_CXX_FLAGS_RELEASE}
+ -DCMAKE_CXX_FLAGS_RELWITHDEBINFO:STRING=${CMAKE_CXX_FLAGS_RELWITHDEBINFO}
+ -DCMAKE_INSTALL_PREFIX:PATH=${INSTALL_PREFIX}
+ -DCMAKE_BUILD_TYPE:STRING=${build_type}
+)
+SET(GDCM_DIR ${build_prefix}/GDCM)
+#=========================================================
+
#=========================================================
# ITK
ExternalProject_Add(
ITK
SOURCE_DIR ${source_prefix}/itk
GIT_REPOSITORY git://itk.org/ITK.git
- GIT_TAG v4.2.0
+ GIT_TAG v4.4.0
INSTALL_COMMAND ""
CMAKE_ARGS
-DCMAKE_C_FLAGS_DEBUG:STRING=${CMAKE_C_FLAGS_DEBUG}
ExternalProject_Add(
VV
- DEPENDS QT VTK ITK
+ DEPENDS QT VTK ITK GDCM
SOURCE_DIR ${source_prefix}/vv
GIT_REPOSITORY git://git.creatis.insa-lyon.fr/clitk
- INSTALL_COMMAND ${MAKE_COMMAND} package
+ INSTALL_DIR ${install_prefix}
+ INSTALL_COMMAND make install
CMAKE_ARGS
-DQT_QMAKE_EXECUTABLE:FILEPATH=${qmake_executable}
+ -DGDCM_DIR:PATH=${GDCM_DIR}
-DITK_DIR:PATH=${ITK_DIR}
-DVTK_DIR:PATH=${VTK_DIR}
-DCMAKE_C_FLAGS_DEBUG:STRING=${CMAKE_C_FLAGS_DEBUG}
WRAP_GGO(clitkMedianImageFilter_GGO_C clitkMedianImageFilter.ggo)
ADD_LIBRARY(clitkMedianImageFilterLib clitkMedianImageGenericFilter.cxx ${clitkMedianImageFilter_GGO_C})
+#=========================================================
IF (CLITK_BUILD_TOOLS)
+
+ #=========================================================
+ # install scripts
+ ADD_SUBDIRECTORY(${CLITK_SOURCE_DIR}/cluster_tools ${PROJECT_BINARY_DIR}/cluster_tools)
+
WRAP_GGO(clitkDicomInfo_GGO_C clitkDicomInfo.ggo)
ADD_EXECUTABLE(clitkDicomInfo clitkDicomInfo.cxx ${clitkDicomInfo_GGO_C})
TARGET_LINK_LIBRARIES(clitkDicomInfo clitkCommon)
WRAP_GGO(clitkMedianTemporalDimension_GGO_C clitkMedianTemporalDimension.ggo)
ADD_EXECUTABLE(clitkMedianTemporalDimension clitkMedianTemporalDimension.cxx
- ${clitkMedianTemporalDimension_GGO_C})
+ ${clitkMedianTemporalDimension_GGO_C})
TARGET_LINK_LIBRARIES(clitkMedianTemporalDimension clitkCommon )
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkMedianTemporalDimension)
TARGET_LINK_LIBRARIES(clitkElastixTransformToMatrix clitkCommon )
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkElastixTransformToMatrix)
+ WRAP_GGO(clitkMatrixToElastixTransform_GGO_C clitkMatrixToElastixTransform.ggo)
+ ADD_EXECUTABLE(clitkMatrixToElastixTransform clitkMatrixToElastixTransform.cxx ${clitkMatrixToElastixTransform_GGO_C})
+ TARGET_LINK_LIBRARIES(clitkMatrixToElastixTransform clitkCommon )
+ SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkMatrixToElastixTransform)
+
+ WRAP_GGO(clitkMatrixInverse_GGO_C clitkMatrixInverse.ggo)
+ ADD_EXECUTABLE(clitkMatrixInverse clitkMatrixInverse.cxx ${clitkMatrixInverse_GGO_C})
+ TARGET_LINK_LIBRARIES(clitkMatrixInverse clitkCommon )
+ SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkMatrixInverse)
WRAP_GGO(clitkSetBackground_GGO_C clitkSetBackground.ggo)
ADD_EXECUTABLE(clitkSetBackground clitkSetBackground.cxx clitkSetBackgroundGenericFilter.cxx ${clitkSetBackground_GGO_C})
TARGET_LINK_LIBRARIES(clitkSetBackground clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkSetBackground)
- WRAP_GGO(clitkGammaIndex_GGO_C clitkGammaIndex.ggo)
+ WRAP_GGO(clitkGammaIndex_GGO_C clitkGammaIndex.ggo)
ADD_EXECUTABLE(clitkGammaIndex clitkGammaIndex.cxx ${clitkGammaIndex_GGO_C})
TARGET_LINK_LIBRARIES(clitkGammaIndex clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkGammaIndex)
WRAP_GGO(clitkUnsharpMask_GGO_C clitkUnsharpMask.ggo)
ADD_EXECUTABLE(clitkUnsharpMask clitkUnsharpMask.cxx ${clitkUnsharpMask_GGO_C})
- TARGET_LINK_LIBRARIES(clitkUnsharpMask clitkCommon )
+ TARGET_LINK_LIBRARIES(clitkUnsharpMask clitkCommon )
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkUnsharpMask)
WRAP_GGO(clitkFooImage_GGO_C clitkFooImage.ggo)
ADD_EXECUTABLE(clitkFooImage clitkFooImage.cxx ${clitkFooImage_GGO_C})
- TARGET_LINK_LIBRARIES(clitkFooImage clitkCommon )
+ TARGET_LINK_LIBRARIES(clitkFooImage clitkCommon )
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkFooImage)
#WRAP_GGO(clitkMedianImageFilter_GGO_C clitkMedianImageFilter.ggo)
# ADD_EXECUTABLE(clitkExtractSlice clitkExtractSlice.cxx clitkExtractSliceGenericFilter.cxx ${clitkExtractSlice_GGO_C})
# TARGET_LINK_LIBRARIES(clitkExtractSlice clitkCommon)
#SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkExtractSlice)
-
+
WRAP_GGO(clitkFlipImage_GGO_C clitkFlipImage.ggo)
ADD_EXECUTABLE(clitkFlipImage clitkFlipImage.cxx clitkFlipImageGenericFilter.cxx ${clitkFlipImage_GGO_C})
TARGET_LINK_LIBRARIES(clitkFlipImage clitkCommon)
TARGET_LINK_LIBRARIES(clitkTransformLandmarks clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkTransformLandmarks)
+ WRAP_GGO(clitkDice_GGO_C clitkDice.ggo)
+ ADD_EXECUTABLE(clitkDice clitkDice.cxx ${clitkDice_GGO_C})
+ TARGET_LINK_LIBRARIES(clitkDice clitkSegmentationGgoLib clitkCommon ${ITK_LIBRARIES} )
+ SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDice)
+
WRAP_GGO(clitkMaskLandmarks_GGO_C clitkMaskLandmarks.ggo)
ADD_EXECUTABLE(clitkMaskLandmarks clitkMaskLandmarks.cxx ${clitkMaskLandmarks_GGO_C})
TARGET_LINK_LIBRARIES(clitkMaskLandmarks clitkCommon)
ADD_EXECUTABLE(clitkJacobianImage clitkJacobianImage.cxx clitkJacobianImageGenericFilter.cxx ${clitkJacobianImage_GGO_C})
TARGET_LINK_LIBRARIES(clitkJacobianImage clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkJacobianImage)
-
+
WRAP_GGO(clitkPadImage_GGO_C clitkPadImage.ggo)
ADD_EXECUTABLE(clitkPadImage clitkPadImage.cxx clitkPadImageGenericFilter.cxx ${clitkPadImage_GGO_C})
TARGET_LINK_LIBRARIES(clitkPadImage clitkCommon)
WRAP_GGO(clitkAnisotropicDiffusion_GGO_C clitkAnisotropicDiffusion.ggo)
ADD_EXECUTABLE(clitkAnisotropicDiffusion clitkAnisotropicDiffusion.cxx
- clitkAnisotropicDiffusionGenericFilter.cxx
- ${clitkAnisotropicDiffusion_GGO_C})
+ clitkAnisotropicDiffusionGenericFilter.cxx
+ ${clitkAnisotropicDiffusion_GGO_C})
TARGET_LINK_LIBRARIES(clitkAnisotropicDiffusion clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkAnisotropicDiffusion)
WRAP_GGO(clitkChangeImageInfo_GGO_C clitkChangeImageInfo.ggo)
ADD_EXECUTABLE(clitkChangeImageInfo clitkChangeImageInfo.cxx
- clitkChangeImageInfoGenericFilter.cxx
- ${clitkChangeImageInfo_GGO_C})
+ clitkChangeImageInfoGenericFilter.cxx
+ ${clitkChangeImageInfo_GGO_C})
TARGET_LINK_LIBRARIES(clitkChangeImageInfo clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkChangeImageInfo)
+ WRAP_GGO(clitkMergeAsciiDoseActor_GGO_C clitkMergeAsciiDoseActor.ggo)
+ ADD_EXECUTABLE(clitkMergeAsciiDoseActor clitkMergeAsciiDoseActor.cxx ${clitkMergeAsciiDoseActor_GGO_C})
+ TARGET_LINK_LIBRARIES(clitkMergeAsciiDoseActor ITKCommon clitkCommon)
+ SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkMergeAsciiDoseActor)
+
+ WRAP_GGO(clitkImageUncertainty_GGO_C clitkImageUncertainty.ggo)
+ ADD_EXECUTABLE(clitkImageUncertainty clitkImageUncertainty.cxx clitkImageUncertainty_ggo.c)
+ TARGET_LINK_LIBRARIES(clitkImageUncertainty clitkCommon ${ITK_LIBRARIES})
+ SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkImageUncertainty)
+
+ #=========================================================
+ option(CLITK_USE_ROOT "Build experimental tools using root" OFF)
+ if (CLITK_USE_ROOT)
+ FIND_PACKAGE(ROOT REQUIRED)
+ if(ROOT_FOUND)
+ MESSAGE(STATUS "ROOT found : ${ROOT_LIBRARY_DIR} ${ROOT_INCLUDE_DIR} ${ROOT_LIBRARIES}")
+ ELSE(ROOT_FOUND)
+ MESSAGE(FATAL_ERROR
+ "Cannot build without ROOT. Please set ROOTSYS environement variable.")
+ endif(ROOT_FOUND)
+ INCLUDE_DIRECTORIES(${ROOT_INCLUDE_DIR})
+ LINK_DIRECTORIES(${ROOT_LIBRARY_DIR})
+ WRAP_GGO(clitkMergeRootFiles_GGO_C clitkMergeRootFiles.ggo)
+ ADD_EXECUTABLE(clitkMergeRootFiles clitkMergeRootFiles.cxx GateMergeManager.cc ${clitkMergeRootFiles_GGO_C})
+ TARGET_LINK_LIBRARIES(clitkMergeRootFiles ${ROOT_LIBRARIES})
+ SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkMergeRootFiles)
+ endif()
+ #=========================================================
+
+
+ #=========================================================
IF(CLITK_EXPERIMENTAL)
WRAP_GGO(clitkBinaryImageToMesh_GGO_C clitkBinaryImageToMesh.ggo)
ADD_EXECUTABLE(clitkBinaryImageToMesh clitkBinaryImageToMesh.cxx ${clitkBinaryImageToMesh_GGO_C})
TARGET_LINK_LIBRARIES(clitkMeshViewer clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkMeshViewer)
ENDIF(CLITK_EXPERIMENTAL)
+ #=========================================================
+
+ #=========================================================
IF(ITK_VERSION_MAJOR VERSION_LESS 4)
MESSAGE("clitkDicomRTPlan2Gate is not compatible with GDCM<2 (ITK<4). It will not be built.")
ELSE(ITK_VERSION_MAJOR VERSION_LESS 4)
TARGET_LINK_LIBRARIES(clitkDicomRTPlan2Gate clitkCommon)
SET(TOOLS_INSTALL ${TOOLS_INSTALL} clitkDicomRTPlan2Gate)
ENDIF(ITK_VERSION_MAJOR VERSION_LESS 4)
+ #=========================================================
+
INSTALL (TARGETS ${TOOLS_INSTALL} DESTINATION bin PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE)
ENDIF(CLITK_BUILD_TOOLS)
-
--- /dev/null
+/*----------------------
+ GATE version name: gate_v...
+
+ Copyright (C): OpenGATE Collaboration
+
+This software is distributed under the terms
+of the GNU Lesser General Public Licence (LGPL)
+See GATE/LICENSE.txt for further details
+----------------------*/
+
+
+#include <TROOT.h>
+#include <TFile.h>
+#include <TChain.h>
+#include <TTree.h>
+#include <TBranch.h>
+#include <TIterator.h>
+#include <TObject.h>
+#include <TKey.h>
+#include <TH1.h>
+#include <TH2.h>
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <glob.h>
+#include <ctype.h>
+#include <sstream>
+#include <vector>
+#include <cstdlib>
+#include <cmath>
+#include <list>
+
+#include "GateMergeManager.hh"
+
+using namespace std;
+
+void GateMergeManager::StartMerging(string splitfileName){
+
+ // get the files to merge
+ ReadSplitFile(splitfileName);
+ //do the merging
+ if (m_fastMerge==true) FastMergeRoot();
+ else MergeRoot();
+
+ //if we are here the merging has been successful
+ //we mark the directory as ready for cleanup
+ //string dir = splitfileName.substr(0,splitfileName.rfind("/"));
+ //string ready = "touch "+dir+"/ready_for_delete";
+ //const int res=system(ready.c_str());
+ //if(res) cout<<"Strange?? Can't mark "<<dir<<" as done!"<<endl;
+};
+
+void GateMergeManager::StartMergingFromFilenames(Strings filenames, string outputfile)
+{
+ for (Strings::const_iterator iter=filenames.begin(); iter!=filenames.end(); iter++)
+ {
+ m_vRootFileNames.push_back(*iter);
+ if(m_verboseLevel>2) cout<<"Root input file name: "<<m_vRootFileNames.back()<<endl;
+ }
+
+ m_Nfiles = m_vRootFileNames.size();
+ m_RootTargetName=m_outDir+outputfile;
+ if(m_verboseLevel>2) cout<<"Root output file name: "<<m_RootTargetName<<endl;
+
+ if (m_fastMerge==true) FastMergeRoot();
+ else MergeRoot();
+};
+
+// to process the splitfile
+void GateMergeManager::ReadSplitFile(string splitfileName){
+
+ ifstream splitfile;
+ splitfile.open(splitfileName.c_str());
+ if(!splitfile){
+ cout<<"Can't open split file: "<<splitfileName<<endl;
+ exit(0);
+ };
+
+ //avoid multiple calls
+ if(m_vRootFileNames.size()>0) return;
+
+ stringstream ssnfiles;
+ char cline[512];
+ while(splitfile){ //we allow for rubbish before "Number of files:"
+ splitfile.getline(cline,512);
+ if(!strncmp(cline,"Number of files:",16)){
+ ssnfiles.str(cline+16);
+ ssnfiles>>m_Nfiles; // the number of files to merge
+ break;
+ }
+ }
+
+ // now we look for the file names
+ int iRoot=0;
+ while(splitfile){
+ splitfile.getline(cline,512);
+ // check for root
+ // input files
+ if(!strncmp(cline,"Root filename:",14)){
+ m_vRootFileNames.push_back(strtok(cline+15," "));
+ m_vRootFileNames[iRoot]+=".root";
+ if(m_verboseLevel>2) cout<<"Root input file name: "<<m_vRootFileNames[iRoot]<<endl;
+ iRoot++;
+ }
+ // output file
+ else if(!strncmp(cline,"Original Root filename:",23)){
+ m_RootTargetName=strtok(cline+23," ");
+ m_RootTargetName+=".root";
+ if(m_outDir!=""){// output directory option used
+ //replace path with outDir
+ size_t pos=m_RootTargetName.rfind('/',m_RootTargetName.length());
+ if(pos==string::npos) pos=-1;
+ m_RootTargetName=m_outDir+m_RootTargetName.substr(pos+1);
+ }
+ if(m_verboseLevel>2) cout<<"Root output file name: "<<m_RootTargetName<<endl;
+ }
+ }
+
+ // check if number of root files correct
+ if(iRoot && iRoot!=m_Nfiles) {
+ cout<<"Inconsistent number of root file entries in split file!"<<endl;
+ exit(0);
+ }
+}
+
+/************************************************************************************/
+void GateMergeManager::FastMergeRoot()
+{
+ //check target name
+ if(m_RootTargetName=="") {
+ cout<<"No root output file name given in split file"<<endl;
+ exit(0);
+ }
+
+ // number of files to merge
+ int nfiles=m_vRootFileNames.size();
+ if(nfiles<1){
+ cout<<"No files to merge!"<<endl;
+ exit(0);
+ }
+ //we try to recover the last_event_ID in all root files
+ for(int i=0;i<nfiles;i++) {
+ filearr.push_back(TFile::Open(m_vRootFileNames[i].c_str(),"OLD"));
+ if(filearr[i]==NULL){
+ cout<<"Not a readable file "<<m_vRootFileNames[i]<<" - exit!"<<endl;
+ exit(0);
+ }
+ }
+ m_lastEvents.resize(nfiles+1);
+ m_lastEvents[0]=0;
+ bool flgLstEvntID(false);
+
+ // latest_event_ID histogram
+ flgLstEvntID=true;
+ TH1* h=NULL;
+ for(int i=0;i<nfiles;i++){
+ h = (TH1*)filearr[i]->Get("latest_event_ID");
+ m_lastEvents[i+1]=int(h->GetMean());
+ }
+
+ if(!flgLstEvntID) {
+ cout<<"FastMerge: No latest_event_ID histogram - exit!"<<endl;
+ exit(0);
+ }
+
+ /*for (int i=0;i<nfiles;i++)
+ {
+ cout<<"Last ID is: "<<m_lastEvents[i+1]<<" File " <<i+1<<endl;
+ }*/
+
+ //find all trees
+ vector<string> treeNames;
+ TFile* node = TFile::Open(m_vRootFileNames[0].c_str(),"OLD");
+ TIter nextkey(node->GetListOfKeys());
+ TKey *key=0;
+ while ((key = (TKey*)nextkey()))
+ {
+ node->cd();
+ TObject* obj = (TObject*)key->ReadObj();
+ if( obj->IsA()->InheritsFrom("TTree"))
+ {
+ //no doubles
+ bool singleName=true;
+ for(unsigned int i=0;i<treeNames.size();i++)
+ {
+ if(treeNames[i].compare(obj->GetName())==0)
+ {
+ singleName=false;
+ break;
+ }
+ }
+ if(singleName) treeNames.push_back(obj->GetName());
+ }
+ }
+ //take care of all trees
+ for(unsigned int i=0;i<treeNames.size();i++)
+ {
+ if(!MergeTree(treeNames[i])) if(m_verboseLevel>1) cout<<"Problem with merging "<<treeNames[i]<<endl;
+ }
+}
+/************************************************************************************/
+void GateMergeManager::MergeRoot(){
+
+ //check target name
+ if(m_RootTargetName=="") {
+ cout<<"No root output file name given in split file"<<endl;
+ exit(0);
+ }
+
+ // number of files to merge
+ int nfiles=m_vRootFileNames.size();
+
+ // the root output file
+ if(m_forced) {
+ m_RootTarget = TFile::Open(m_RootTargetName.c_str(),"RECREATE");
+ } else {
+ m_RootTarget = TFile::Open(m_RootTargetName.c_str(),"NEW");
+ if(!m_RootTarget){
+ cout<<"The root ouput file already exists! Try -f to overwrite and erase the file."<<endl;
+ exit(0);
+ }
+ }
+
+ if(nfiles<1){
+ cout<<"No files to merge!"<<endl;
+ exit(0);
+ }
+
+ TFile* filearr[nfiles];
+ for(int i=1;i<nfiles;i++) {
+ filearr[i] = TFile::Open(m_vRootFileNames[i].c_str(),"OLD");
+ if(filearr[i]==NULL){
+ cout<<"Not a readable file "<<m_vRootFileNames[i]<<" - exit!"<<endl;
+ exit(0);
+ }
+ }
+
+// first we copy all histos (only top directory)
+// and look for the latestEventID
+// and find out which trees/ntuples exist
+ m_lastEvents.resize(nfiles+1);
+ m_lastEvents[0]=0;
+
+ bool flgLstEvntID(false);
+
+ if(m_verboseLevel>0) {
+ cout<<"Combining: ";
+ for(int i=0;i<nfiles;i++) cout<<m_vRootFileNames[i]<<" ";
+ cout<<"-> "<<m_RootTarget->GetName()<<endl;
+ }
+
+ vector<string> treeNames;
+ TFile* node = TFile::Open(m_vRootFileNames[0].c_str(),"OLD");
+ if(node==NULL){
+ cout<<"Not a readable file "<<m_vRootFileNames[0]<<" - exit!"<<endl;
+ exit(0);
+ }
+ TIter nextkey(node->GetListOfKeys());
+ TKey *key=0;
+ while ((key = (TKey*)nextkey())) {
+ node->cd();
+ TObject* obj = (TObject*)key->ReadObj();
+ if( obj->IsA()->InheritsFrom("TH1")){
+ m_RootTarget->cd();
+ TH1 *h1 = (TH1 *)obj->Clone();
+ if(strcmp(h1->GetName(),"latest_event_ID")){
+ //any other histogram
+ for(int i=1;i<nfiles;i++){
+ TH1 *h2 = (TH1*)filearr[i]->Get( h1->GetName() );
+ h1->Add( h2 );
+ }
+ h1->Write();
+ } else {
+ // latest_event_ID histogram
+ flgLstEvntID=true;
+ m_lastEvents[1]=int(h1->GetMean());
+ TH1* h=NULL;
+ for(int i=1;i<nfiles;i++){
+ h = (TH1*)filearr[i]->Get("latest_event_ID");
+ m_lastEvents[i+1]=int(h->GetMean());
+ }
+ if(h!=NULL) h->Write(); // we keep the highest latest_event_ID to allow for successive merging
+ else h1->Write();
+ }
+ }
+ if( obj->IsA()->InheritsFrom("TH2")){
+ m_RootTarget->cd();
+ TH2 *h1 = (TH2 *)obj->Clone();
+ for(int i=1;i<nfiles;i++){
+ TH2 *h2 = (TH2*)filearr[i]->Get( h1->GetName() );
+ h1->Add( h2 );
+ }
+ h1->Write();
+ }
+ if( obj->IsA()->InheritsFrom("TTree")){
+ //no doubles
+ bool singleName=true;
+ for(unsigned int i=0;i<treeNames.size();i++) {
+ if(treeNames[i].compare(obj->GetName())==0){
+ singleName=false;
+ break;
+ }
+ }
+ if(singleName) treeNames.push_back(obj->GetName());
+ }
+ }
+ if(!flgLstEvntID) {
+ cout<<"Merge: No latest_event_ID histogram - exit!"<<endl;
+ exit(0);
+ }
+
+ //now we take care of the trees
+ for(unsigned int i=0;i<treeNames.size();i++)
+ if(!MergeTree(treeNames[i])) if(m_verboseLevel>1) cout<<"Problem with merging "<<treeNames[i]<<endl;
+
+ // everything is done
+};
+
+/*******************************************************************************************/
+// cleanup after merging
+void GateMergeManager::StartCleaning(string splitfileName,bool test){
+
+ // get the files to erase
+ ReadSplitFile(splitfileName);
+
+ string dir=splitfileName.substr(0,splitfileName.rfind("/"));
+
+ // test for the mark: ready_for_delete
+ string touched=dir+"/ready_for_delete";
+ ifstream ready;
+ ready.open(touched.c_str());
+ if (!ready) {
+ cout<<"Cannot do the cleanup - directory "<<dir<<" not marked as ready!"<<endl;
+ if(!test) exit(0);
+ }
+
+ // print info
+ if(test) cout<<"I would like to erase:"<<endl;
+ if(!test&&m_verboseLevel>1) cout<<"Going to erase the following files:"<<endl;
+ if(m_verboseLevel>1||test){
+ cout<<" --> "<<dir<<endl;
+ // root
+ for(unsigned int i=0;i<m_vRootFileNames.size();i++){
+ cout<<" --> "<<m_vRootFileNames[i]<<endl;
+ }
+ }
+
+ // erase!
+ if(!test){
+ for(unsigned int i=0;i<m_vRootFileNames.size();i++){
+ const string rmfiles("rm -f "+m_vRootFileNames[i]);
+ const int res = system(rmfiles.c_str());
+ if(res) {
+ cout<<"Could not remove "<<m_vRootFileNames[i]<<endl;
+ cout<<"Please remove it manually!"<<endl;
+ }
+ }
+ const string rmfiles("rm -f "+dir+"/*");
+ system(rmfiles.c_str());
+ const string rmdir="rm -f -r "+dir;
+ const int res2 = system(rmdir.c_str());
+ if(res2) {
+ cout<<"Could not remove "<<dir<<endl;
+ cout<<"Please remove it manually!"<<endl;
+ }
+ }
+}
+
+/*******************************************************************************************/
+//find out what kind of tree we have and call the right merger
+bool GateMergeManager::MergeTree(string chainName){
+if (m_fastMerge==false)
+ {
+ TChain* chain = new TChain(chainName.c_str());
+
+ // number of files to merge
+ int nfiles=m_vRootFileNames.size();
+
+ for(int i=0;i<nfiles;i++) chain->Add(m_vRootFileNames[i].c_str());
+ int nentries=chain->GetEntries();
+ if(nentries<=0) {
+ if(m_verboseLevel>1) cout<<chain->GetName()<<" is empty!"<<endl;
+ return false;
+ }
+
+ if(chainName=="Gate") MergeGate(chain);
+
+ if(chain->FindBranch("eventID1")!=NULL) {
+ if( (chain->FindBranch("runID")==NULL)
+ ||(chain->FindBranch("time1")==NULL)
+ ||(chain->FindBranch("time2")==NULL) ) {
+ cout<<"Cannot find one of: runID, time1, time2 in "<<chain->GetName()<<endl;
+ return false;
+ }
+ // coincidence
+ return MergeCoin(chain);
+
+ } else if( (chain->FindBranch("runID")==NULL)||(chain->FindBranch("eventID")==NULL)
+ ||(chain->FindBranch("time")==NULL) )
+ {
+ cout<<"Cannot find one of: runID, eventID, time in "<<chain->GetName()<<endl;
+ return false;
+ }
+ // Singles or Hits
+ return MergeSing(chain);
+ }
+ else
+ {
+ if (chainName.find("Hits",0)!=string::npos) return FastMergeSing(chainName);
+ if (chainName.find("Singles",0)!=string::npos) return FastMergeSing(chainName);
+ if (chainName.find("Gate",0)!=string::npos) return FastMergeGate(chainName);
+ if (chainName.find("Coincidences",0)!=string::npos) return FastMergeCoin(chainName);
+ }
+ return 0;
+}
+
+/*******************************************************************************************/
+bool GateMergeManager::FastMergeGate(string name)
+{
+//create the output file
+float event = 0;
+int offset = 0;
+int currentfile = 0;
+float lastEvent =-1;
+//float maxtime =-999999999;
+string clusterName=name+"_cluster";
+
+for(int j=0;j<m_Nfiles;j++)
+ {
+ filearr[j]->ReOpen("UPDATE");
+ TTree *oldTree = (TTree*)gDirectory->Get(name.c_str());
+ TBranch *branch = oldTree->GetBranch("event");
+ branch->SetAddress(&event);
+ TBranch* newBranch=oldTree->Branch("eventcluster",&event,"event");
+
+ Int_t nentries=(Int_t)oldTree->GetEntries();
+
+ if (j==0) offset=0;
+ else offset+=m_lastEvents[currentfile];
+ currentfile++;
+
+ for(int i=0;i<nentries;i++)
+ {
+ branch->GetEvent(i);
+ if(lastEvent!=event)
+ {
+ lastEvent=event;
+ event+=offset;
+ newBranch->Fill();
+ }
+ else
+ {
+ if((i==nentries-1) && (j==m_Nfiles-1))
+ {
+ event+=offset;
+ newBranch->Fill();
+ }
+ else if (i!=nentries-1)
+ {
+ branch->GetEvent(i);
+ event+=offset;
+ newBranch->Fill();
+ lastEvent=0;
+ offset=0;
+ }
+ }
+ }
+ oldTree->Write();
+ }
+ return true;
+}
+
+/*******************************************************************************************/
+bool GateMergeManager::FastMergeSing(string name)
+{
+int eventID = 0;
+int runID = 0;
+int offset = 0;
+int currentfile = 0;
+float lastRun =-1;
+string clusterName=name+"_cluster";
+
+for(int j=0;j<m_Nfiles;j++)
+ {
+ filearr[j]->ReOpen("UPDATE");
+ TTree *oldTree = (TTree*)gDirectory->Get(name.c_str());
+ TBranch *branch = oldTree->GetBranch("eventID");
+ branch->SetAddress(&eventID);
+ TBranch *branch2 = oldTree->GetBranch("runID");
+ branch2->SetAddress(&runID);
+ TBranch* newBranch=oldTree->Branch("eventIDcluster",&eventID,"eventID/I");
+
+ Int_t nentries=(Int_t)oldTree->GetEntries();
+
+ if (j==0) offset=0;
+ else offset+=m_lastEvents[currentfile];
+ currentfile++;
+
+ for(int i=0;i<nentries;i++)
+ {
+ branch->GetEvent(i);
+ if(lastRun!=runID)
+ {
+ lastRun=runID;
+ offset=0;
+ }
+ eventID+=offset;
+ newBranch->Fill();
+ }
+ oldTree->Write();
+ }
+ return true;
+}
+
+/*******************************************************************************************/
+bool GateMergeManager::FastMergeCoin(string name)
+{
+//create the output file
+int eventID1 = 0;
+int eventID2 = 0;
+int runID = 0;
+int offset = 0;
+int currentfile = 0;
+float lastRun =-1;
+string clusterName=name+"_cluster";
+cout<<"starting coincidence merging..."<<endl;
+for(int j=0;j<m_Nfiles;j++)
+ {
+cout<<"working on file..."<<j<<endl;
+ filearr[j]->ReOpen("UPDATE");
+ TTree *oldTree = (TTree*)gDirectory->Get(name.c_str());
+ TBranch *branch1 = oldTree->GetBranch("eventID1");
+ branch1->SetAddress(&eventID1);
+ TBranch *branch2 = oldTree->GetBranch("eventID2");
+ branch2->SetAddress(&eventID2);
+ TBranch *branch3 = oldTree->GetBranch("runID");
+ branch3->SetAddress(&runID);
+
+ TBranch* newBranch1=oldTree->Branch("eventID1cluster",&eventID1,"eventID1/I");
+ TBranch* newBranch2=oldTree->Branch("eventID2cluster",&eventID2,"eventID2/I");
+
+ Int_t nentries=(Int_t)oldTree->GetEntries();
+
+ if (j==0) offset=0;
+ else offset+=m_lastEvents[currentfile];
+ cout<<"the offset in this file is: "<<offset<<endl;
+ currentfile++;
+
+ for(int i=0;i<nentries;i++){
+ branch1->GetEvent(i);
+ branch2->GetEvent(i);
+ branch3->GetEvent(i);
+ if(lastRun!=runID)
+ {
+ lastRun=runID;
+ offset=0;
+ }
+ eventID1+=offset;
+ eventID2+=offset;
+ newBranch1->Fill();
+ newBranch2->Fill();
+ }
+ oldTree->Write();
+ }
+ return true;
+}
+
+/*******************************************************************************************/
+// Gate tree merger
+bool GateMergeManager::MergeGate(TChain* chainG) {
+
+ int nentries=chainG->GetEntries();
+
+ float event = 0;
+ chainG->SetBranchAddress("event",&event);
+ Float_t iontime = 0;
+ chainG->SetBranchAddress("iontime",&iontime);
+
+ //create the new tree
+ //Allow for large files and do not write every bit separately
+ m_RootTarget->cd();
+ TTree * newTree = chainG->CloneTree(0);
+ newTree->SetAutoSave(2000000000);
+ if(m_maxRoot!=0) newTree->SetMaxTreeSize(m_maxRoot);
+ else newTree->SetMaxTreeSize(17179869184LL);
+
+ // changing the compression level everywhere
+ TBranch *br;
+ TIter next(newTree->GetListOfBranches());
+ while ((br=(TBranch*)next())) br->SetCompressionLevel(m_CompLevel);
+
+ // the main part: modify the eventID
+ int offset = 0;
+ int currentTree = 0;
+ float lastEvent =-1;
+ float maxtime =-999999999;
+ for(int i=0;i<nentries;i++){
+ if (chainG->GetEntry(i) <= 0) break; //end of chain
+ if(chainG->GetTreeNumber()>currentTree) {
+ currentTree++;
+ offset+=m_lastEvents[currentTree];
+ // check for overlaping time intervalls between different files
+ // (not within the same file i.e. no time order assumed)
+ if(iontime<maxtime)
+ if(m_verboseLevel>0)
+ cout<<"Warning - overlapping Gate iontime ("
+ <<iontime<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
+ }
+ // run end is marked by repeating event
+ if(lastEvent!=event) {
+ lastEvent=event;
+// if (fpclassify(iontime)!=FP_NAN && fpclassify(iontime)!=FP_INFINITE) {
+ if(finite(iontime)){
+ if(maxtime<iontime) maxtime=iontime;
+ // the offset to get a unique event numbering
+ event+=offset;
+ newTree->Fill();
+ } else {
+ cout<<"Warning - inf or NaN in Gate tree for iontime! "<<m_vRootFileNames[currentTree].c_str()<<endl;
+ //we are lost anyhow
+ maxtime =-999999999;
+ }
+ } else {
+ // run end
+ if(chainG->GetEntry(i+1) <= 0) { //chain end fill the double event
+ event+=offset;
+ newTree->Fill();
+ } else if(chainG->GetTreeNumber()==currentTree //no new file or
+ ||(chainG->GetTreeNumber()>currentTree
+ &&m_lastEvents[chainG->GetTreeNumber()]==0) ) { //new file with no offset fill double event
+ chainG->GetEntry(i);
+ event+=offset;
+ newTree->Fill();
+ lastEvent=0;// we may have triple 1 (1 event run)
+ offset=0;
+ }
+ }
+ }
+ newTree->Write();
+ delete newTree;
+ return true;
+}
+
+/*******************************************************************************************/
+// Singles and Hits tree merger
+bool GateMergeManager::MergeSing(TChain* chainS){
+
+ int nentriesS=chainS->GetEntries();
+
+ int eventID = 0;
+ int runID = 0;
+ double time = 0;
+
+ chainS->SetBranchAddress("eventID",&eventID);
+ chainS->SetBranchAddress("runID",&runID);
+ chainS->SetBranchAddress("time",&time);
+
+ m_RootTarget->cd();
+ TTree * newSing = chainS->CloneTree(0);
+ newSing->SetAutoSave(2000000000);
+ if(m_maxRoot!=0) newSing->SetMaxTreeSize(m_maxRoot);
+ else newSing->SetMaxTreeSize(17179869184LL);
+
+ // changing CompLevel everywhere
+ TBranch *br;
+ TIter next(newSing->GetListOfBranches());
+ while ((br=(TBranch*)next())) br->SetCompressionLevel(m_CompLevel);
+
+ //the main part: changing eventID
+ int offset = 0;
+ int currentTree = 0;
+ float lastRun = 0;
+ float maxtime =-999999999;
+ for(int i=0;i<nentriesS;i++){
+ if(chainS->GetEntry(i) <= 0) break; //end of chain
+ if(chainS->GetTreeNumber()>currentTree) {
+ currentTree++;
+ offset+=m_lastEvents[currentTree];
+ // check for overlaping time intervalls between different files
+ // (not within the same file i.e. no time order assumed)
+ if(time<maxtime)
+ if(m_verboseLevel>0)
+ cout<<"Warning - overlapping Singles time ("
+ <<time<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
+ }
+ // the offset to get a unique event numbering
+ if(lastRun!=runID) {
+ // run end
+ lastRun=runID;
+ offset=0; //new run in file we must not change the eventID anymore
+ }
+ eventID+=offset;
+ newSing->Fill();
+ if(maxtime<time)maxtime=time;
+ }
+ newSing->Write();
+ delete newSing;
+ return true;
+};
+
+/*******************************************************************************************/
+// Coincidences tree merger
+bool GateMergeManager::MergeCoin(TChain* chainC){
+
+ int nentriesC=chainC->GetEntries();
+
+ int eventID1 = 0;
+ chainC->SetBranchAddress("eventID1",&eventID1);
+ double time1 = 0;
+ chainC->SetBranchAddress("time1",&time1);
+ int eventID2 = 0;
+ chainC->SetBranchAddress("eventID2",&eventID2);
+ double time2 = 0;
+ chainC->SetBranchAddress("time2",&time2);
+ int runID = 0;
+ chainC->SetBranchAddress("runID",&runID);
+
+ m_RootTarget->cd();
+ TTree * newCoin = chainC->CloneTree(0);
+ newCoin->SetAutoSave(2000000000);
+ if(m_maxRoot!=0) newCoin->SetMaxTreeSize(m_maxRoot);
+ else newCoin->SetMaxTreeSize(17179869184LL);
+
+ // changing CompLevel everywhere
+ TBranch *br;
+ TIter next(newCoin->GetListOfBranches());
+ while ((br=(TBranch*)next())) br->SetCompressionLevel(m_CompLevel);
+
+ int offset = 0;
+ int currentTree = 0;
+ float lastRun = 0;
+ float maxtime =-999999999;
+ for(int i=0;i<nentriesC;i++){
+ if(chainC->GetEntry(i) <= 0) break; //end of chain
+ if(chainC->GetTreeNumber()>currentTree) {
+ currentTree++;
+ offset+=m_lastEvents[currentTree];
+ // check for overlaping time intervalls between different files
+ // (not within the same file i.e. no time order assumed)
+ if(time1<maxtime)
+ if(m_verboseLevel>0)
+ cout<<"Warning - overlapping Coincidences time1 ("
+ <<time1<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
+ if(time2<maxtime)
+ if(m_verboseLevel>0)
+ cout<<"Warning - overlapping Coincidences time2 ("
+ <<time2<<") in file: "<<m_vRootFileNames[currentTree].c_str()<<endl;
+ }
+ // the offset to get a unique event numbering
+ if(lastRun!=runID) {
+ // run end
+ lastRun=runID;
+ offset=0; //new run in file we must not change the eventID anymore
+ }
+ eventID1+=offset;
+ eventID2+=offset;
+ newCoin->Fill();
+ if(maxtime<time1)maxtime=time1;
+ if(maxtime<time2)maxtime=time2;
+ }
+ newCoin->Write();
+ delete newCoin;
+ return true;
+}
+/*******************************************************************************************/
--- /dev/null
+/*----------------------
+ GATE version name: gate_v...
+
+ Copyright (C): OpenGATE Collaboration
+
+This software is distributed under the terms
+of the GNU Lesser General Public Licence (LGPL)
+See GATE/LICENSE.txt for further details
+----------------------*/
+
+
+#ifndef GateMergeManager_h
+#define GateMergeManager_h 1
+#include <iostream>
+#include <string>
+#include <vector>
+#include <TFile.h>
+#include <TChain.h>
+#include <cstdlib>
+#include <list>
+
+using namespace std;
+typedef list<string> Strings;
+
+class GateMergeManager
+{
+public:
+
+ GateMergeManager(bool fastMerge,int verboseLevel,bool forced,Long64_t maxRoot,string outDir){
+ m_verboseLevel = verboseLevel;
+ m_forced = forced;
+ m_maxRoot = maxRoot;
+ m_outDir = outDir;
+ m_CompLevel = 1;
+ m_fastMerge = fastMerge;
+ };
+ ~GateMergeManager()
+ {
+ for (std::vector<TFile*>::const_iterator iter=filearr.begin(); iter!=filearr.end(); iter++)
+ {
+ (*iter)->Close();
+ delete *iter;
+ }
+ }
+
+
+ void StartMerging(string splitfileName);
+ void StartMergingFromFilenames(Strings filenames, string outputfile);
+ void ReadSplitFile(string splitfileName);
+ bool MergeTree(string name);
+ bool MergeGate(TChain* chain);
+ bool MergeSing(TChain* chain);
+ bool MergeCoin(TChain* chain);
+
+ // the cleanup after succesful merging
+ void StartCleaning(string splitfileName,bool test);
+
+ // the merging methods
+ void MergeRoot();
+
+private:
+ void FastMergeRoot();
+ bool FastMergeGate(string name);
+ bool FastMergeSing(string name);
+ bool FastMergeCoin(string name);
+ bool m_forced; // if to overwrite existing files
+ int m_verboseLevel;
+ std::vector<TFile*> filearr;
+ Long64_t m_maxRoot; // maximum size of root output file
+ int m_CompLevel; // compression level for root output
+ string m_outDir; // where to save the output files
+ int m_Nfiles; // number of files to mergw
+ vector<int> m_lastEvents; // latestevent from al files
+ vector<string> m_vRootFileNames; // names of root files to merge
+ TFile* m_RootTarget; // root output file
+ string m_RootTargetName; // name of target i.e. root output file
+ bool m_fastMerge; // fast merge option, corrects the eventIDs locally
+};
+
+
+#endif
//----------------------------------------
void Update();
- template<unsigned int Dimension, class PixelType>
- static
- typename itk::Matrix<double, Dimension+1, Dimension+1>
- createMatrixFromElastixFile(std::vector<std::string> & filename, bool verbose=true);
-
protected:
//----------------------------------------
template <unsigned int Dimension, class PixelType> void UpdateWithDimAndPixelType();
template <unsigned int Dimension, class PixelType> void UpdateWithDimAndVectorType();
- static bool GetElastixValueFromTag(std::ifstream & is, std::string tag, std::string & value);
- static void GetValuesFromValue(const std::string & s,
- std::vector<std::string> & values);
-
//----------------------------------------
// Data members
//----------------------------------------
#include <istream>
#include <iterator>
#include <itkCenteredEuler3DTransform.h>
+#include "clitkElastix.h"
namespace clitk
{
if (m_ArgsInfo.elastix_given) {
std::vector<std::string> s;
for(uint i=0; i<m_ArgsInfo.elastix_given; i++) s.push_back(m_ArgsInfo.elastix_arg[i]);
- matrix = createMatrixFromElastixFile<Dimension,PixelType>(s, m_Verbose);
+ matrix = createMatrixFromElastixFile<Dimension>(s, m_Verbose);
}
else
matrix.SetIdentity();
}
//-------------------------------------------------------------------
-
-
- //-------------------------------------------------------------------
- template<class args_info_type>
- template<unsigned int Dimension, class PixelType>
- typename itk::Matrix<double, Dimension+1, Dimension+1>
- AffineTransformGenericFilter<args_info_type>::createMatrixFromElastixFile(std::vector<std::string> & filename, bool verbose)
- {
- if (Dimension != 3) {
- FATAL("Only 3D yet" << std::endl);
- }
- typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
-
- itk::CenteredEuler3DTransform<double>::Pointer mat = itk::CenteredEuler3DTransform<double>::New();
- itk::CenteredEuler3DTransform<double>::Pointer previous;
- for(uint j=0; j<filename.size(); j++) {
-
- // Open file
- if (verbose) std::cout << "Read elastix parameters in " << filename[j] << std::endl;
- std::ifstream is;
- clitk::openFileForReading(is, filename[j]);
-
- // Check Transform
- std::string s;
- bool b = GetElastixValueFromTag(is, "Transform ", s);
- if (!b) {
- FATAL("Error must read 'Transform' in " << filename[j] << std::endl);
- }
- if (s != "EulerTransform") {
- FATAL("Sorry only 'EulerTransform'" << std::endl);
- }
-
- // FIXME check
- // (InitialTransformParametersFilename[j] "NoInitialTransform")
-
- // Get CenterOfRotationPoint
- GetElastixValueFromTag(is, "CenterOfRotationPoint ", s); // space is needed
- if (!b) {
- FATAL("Error must read 'CenterOfRotationPoint' in " << filename[j] << std::endl);
- }
- std::vector<std::string> cor;
- GetValuesFromValue(s, cor);
-
- // Get Transformparameters
- GetElastixValueFromTag(is, "TransformParameters ", s); // space is needed
- if (!b) {
- FATAL("Error must read 'TransformParameters' in " << filename[j] << std::endl);
- }
- std::vector<std::string> results;
- GetValuesFromValue(s, results);
-
- // construct a stream from the string
- itk::CenteredEuler3DTransform<double>::ParametersType p;
- p.SetSize(9);
- for(uint i=0; i<3; i++)
- p[i] = atof(results[i].c_str()); // Rotation
- for(uint i=0; i<3; i++)
- p[i+3] = atof(cor[i].c_str()); // Centre of rotation
- for(uint i=0; i<3; i++)
- p[i+6] = atof(results[i+3].c_str()); // Translation
- mat->SetParameters(p);
-
- if (verbose) {
- std::cout << "Rotation (deg) : " << rad2deg(p[0]) << " " << rad2deg(p[1]) << " " << rad2deg(p[2]) << std::endl;
- std::cout << "Center of rot (phy) : " << p[3] << " " << p[4] << " " << p[5] << std::endl;
- std::cout << "Translation (phy) : " << p[6] << " " << p[7] << " " << p[8] << std::endl;
- }
-
- // Compose with previous if needed
- if (j!=0) {
- mat->Compose(previous);
- if (verbose) {
- std::cout << "Composed rotation (deg) : " << rad2deg(mat->GetAngleX()) << " " << rad2deg(mat->GetAngleY()) << " " << rad2deg(mat->GetAngleZ()) << std::endl;
- std::cout << "Composed center of rot (phy) : " << mat->GetCenter() << std::endl;
- std::cout << "Compsoed translation (phy) : " << mat->GetTranslation() << std::endl;
- }
- }
- // previous = mat->Clone(); // ITK4
- previous = itk::CenteredEuler3DTransform<double>::New();
- previous->SetParameters(mat->GetParameters());
- }
-
- mat = previous;
- for(uint i=0; i<3; i++)
- for(uint j=0; j<3; j++)
- matrix[i][j] = mat->GetMatrix()[i][j];
- // Offset is -Rc + t + c
- matrix[0][3] = mat->GetOffset()[0];
- matrix[1][3] = mat->GetOffset()[1];
- matrix[2][3] = mat->GetOffset()[2];
- matrix[3][3] = 1;
-
- return matrix;
- }
-
- //-------------------------------------------------------------------
- template<class args_info_type>
- bool
- AffineTransformGenericFilter<args_info_type>::GetElastixValueFromTag(std::ifstream & is,
- std::string tag,
- std::string & value)
- {
- std::string line;
- is.seekg (0, is.beg);
- while(std::getline(is, line)) {
- unsigned pos = line.find(tag);
- if (pos<line.size()) {
- value=line.substr(pos+tag.size(),line.size()-2);// remove last ')'
- value.erase (std::remove (value.begin(), value.end(), '"'), value.end());
- value.erase (std::remove (value.begin(), value.end(), ')'), value.end());
- return true;
- }
- }
- return false;
- }
- //-------------------------------------------------------------------
-
-
- //-------------------------------------------------------------------
- template<class args_info_type>
- void
- AffineTransformGenericFilter<args_info_type>::GetValuesFromValue(const std::string & s,
- std::vector<std::string> & values)
- {
- std::stringstream strstr(s);
- std::istream_iterator<std::string> it(strstr);
- std::istream_iterator<std::string> end;
- std::vector<std::string> results(it, end);
- values.clear();
- values.resize(results.size());
- for(uint i=0; i<results.size(); i++) values[i] = results[i];
- }
- //-------------------------------------------------------------------
-
} //end clitk
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+
+// clitk
+#include "clitkIO.h"
+#include "clitkDice_ggo.h"
+
+#include "clitkLabelImageOverlapMeasureGenericFilter.h"
+
+//--------------------------------------------------------------------
+int main(int argc, char * argv[]) {
+
+ // Init command line
+ GGO(clitkDice, args_info);
+ CLITK_INIT;
+
+ // Filter
+ typedef clitk::LabelImageOverlapMeasureGenericFilter<args_info_clitkDice> FilterType;
+ FilterType::Pointer filter = FilterType::New();
+
+ filter->SetArgsInfo(args_info);
+
+ try {
+ filter->Update();
+ } catch(std::runtime_error e) {
+ std::cout << e.what() << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ return EXIT_SUCCESS;
+} // This is the end, my friend
+//--------------------------------------------------------------------
--- /dev/null
+#File clitkDice.ggo
+package "clitkDice"
+version "1.0"
+purpose "Compute Dice between labeled images. Background must be 0."
+
+section "General options"
+option "config" - "Config file" string no
+option "verbose" v "Verbose" flag off
+option "imagetypes" - "Display allowed image types" flag off
+option "long" l "Long display (not only dice)" flag off
+
+section "Input"
+option "input1" i "Input mask 1" string yes
+option "input2" j "Input mask 2" string yes
+
+option "label1" - "Label in input1" int no default="1"
+option "label2" - "Label in input2" int no default="1"
// clitk
#include "clitkElastixTransformToMatrix_ggo.h"
#include "clitkAffineTransformGenericFilter.h"
+#include "clitkElastix.h"
+#include "clitkMatrix.h"
//--------------------------------------------------------------------
int main(int argc, char * argv[])
CLITK_INIT;
// Use static fct of AffineTransformGenericFilter
- typedef clitk::AffineTransformGenericFilter<args_info_clitkElastixTransformToMatrix> FilterType;
std::vector<std::string> l;
l.push_back(args_info.input_arg);
- itk::Matrix<double, 4, 4> m =
- FilterType::createMatrixFromElastixFile<3, int>(l, args_info.verbose_flag);
+ itk::Matrix<double, 4, 4> m = clitk::createMatrixFromElastixFile<3>(l, args_info.verbose_flag);
// Print matrix
std::ofstream os;
clitk::openFileForWriting(os, args_info.output_arg);
- for(unsigned int i=0; i<4; i++) {
- for(unsigned int j=0; j<4; j++)
- os << m[i][j] << " ";
- os << std::endl;
- }
- os.close();
+ os << clitk::Get4x4MatrixDoubleAsString(m, 16);
+ os.close();
return EXIT_SUCCESS;
}// end main
--- /dev/null
+/*=========================================================================
+
+ Program: clitk
+ Module: $RCSfile: clitkImageUncertainty.cxx,v $
+ Language: C++
+ Date: $Date: 2011/03/03 15:03:30 $
+ Version: $Revision: 1.3 $
+
+ Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
+ l'Image). All rights reserved. See Doc/License.txt or
+ http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+
+#ifndef CLITKIMAGEUNCERTAINTY_CXX
+#define CLITKIMAGEUNCERTAINTY_CXX
+
+/**
+ =================================================
+ * @file clitkImageUncertainty.cxx
+ * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
+ * @date 04 Jul 2006 14:03:57
+ *
+ * @brief
+ *
+ *
+ =================================================*/
+
+// clitk include
+#include "clitkImageUncertainty_ggo.h"
+#include "clitkImageCommon.h"
+#include "clitkCommon.h"
+#include "clitkPortability.h"
+
+// itk include
+#include "itkImageRegionConstIterator.h"
+#include "itkImageRegionIterator.h"
+
+//====================================================================
+int main(int argc, char * argv[]) {
+
+ // init command line
+ GGO(clitkImageUncertainty, args_info);
+
+ // Declare main types
+ typedef float PixelType;
+ const unsigned int Dimension=3;
+ typedef itk::Image< PixelType, Dimension > ImageType;
+
+ // Read images
+ ImageType::Pointer input =
+ clitk::readImage<ImageType>(args_info.input_arg, args_info.verbose_flag);
+ ImageType::Pointer inputSquared =
+ clitk::readImage<ImageType>(args_info.inputSquared_arg, args_info.verbose_flag);
+
+ // Create Output
+ ImageType::Pointer output = ImageType::New();
+ output->SetRegions(input->GetLargestPossibleRegion());
+ output->CopyInformation(input);
+ output->Allocate();
+
+ // Loop
+ typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
+ ConstIteratorType pi(input, input->GetLargestPossibleRegion());
+ ConstIteratorType pii(inputSquared, inputSquared->GetLargestPossibleRegion());
+ pi.GoToBegin();
+ pii.GoToBegin();
+ typedef itk::ImageRegionIterator<ImageType> IteratorType;
+ IteratorType po(output, output->GetLargestPossibleRegion());
+ po.GoToBegin();
+
+ int NumberOfEvents = args_info.NumberOfEvents_arg;
+ while ( !pi.IsAtEnd() ) {
+ double squared = pii.Get();
+ double mean = pi.Get();
+ double uncert = sqrt((NumberOfEvents*squared - mean*mean) / ((NumberOfEvents-1)*(mean*mean)));
+ if (!IsNormal(uncert)) uncert = 1.;
+ po.Set(uncert);
+ ++pi;
+ ++pii;
+ ++po;
+ }
+// *po = sqrt( (NumberOfEvents*squared - mean*mean) /
+// ((NumberOfEvents-1)*(mean*mean)) );
+// ++po;
+
+
+
+ // Write output image
+ // DD(clitk::GetExtension(args_info.output_arg));
+ if (clitk::GetExtension(args_info.output_arg) != "txt") {
+ clitk::writeImage<ImageType>(output, args_info.output_arg, args_info.verbose_flag);
+ }
+ else {
+ std::ofstream os;
+ clitk::openFileForWriting(os, args_info.output_arg);
+ typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
+ IteratorType pi(output, output->GetLargestPossibleRegion());
+ pi.GoToBegin();
+ os << "# Image size = " << output->GetLargestPossibleRegion().GetSize() << std::endl;
+ os << "# Image spacing = " << output->GetSpacing() << std::endl;
+ os << "# Number of events = " << NumberOfEvents << std::endl;
+ while (!pi.IsAtEnd()) {
+ os << pi.Get() << std::endl;
+ ++pi;
+ }
+ }
+}
+
+
+#endif /* end #define CLITKIMAGEUNCERTAINTY_CXX */
--- /dev/null
+# file clitkImageRescaleIntensity.ggo
+package "clitk"
+version "Rescale intensity in the image"
+
+option "config" - "Config file" string no
+option "input" i "Input image filename" string yes
+option "inputSquared" s "Input squared image filename" string yes
+option "output" o "Output image filename" string yes
+option "NumberOfEvents" n "Number of events" int yes
+option "verbose" v "Verbose" flag off
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+
+#ifndef CLITKLABELIMAGEOVERLAPMEASUREGENERICFILTER_H
+#define CLITKLABELIMAGEOVERLAPMEASUREGENERICFILTER_H
+
+// clitk
+#include "clitkImageToImageGenericFilter.h"
+#include "clitkLabelImageOverlapMeasureFilter.h"
+#include "clitkBoundingBoxUtils.h"
+#include "clitkCropLikeImageFilter.h"
+
+//--------------------------------------------------------------------
+namespace clitk
+{
+
+ template<class ArgsInfoType>
+ class ITK_EXPORT LabelImageOverlapMeasureGenericFilter:
+ public ImageToImageGenericFilter<LabelImageOverlapMeasureGenericFilter<ArgsInfoType> >
+ {
+ public:
+ //--------------------------------------------------------------------
+ LabelImageOverlapMeasureGenericFilter();
+
+ //--------------------------------------------------------------------
+ typedef ImageToImageGenericFilter<LabelImageOverlapMeasureGenericFilter<ArgsInfoType> > Superclass;
+ typedef LabelImageOverlapMeasureGenericFilter Self;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ //--------------------------------------------------------------------
+ itkNewMacro(Self);
+ itkTypeMacro(LabelImageOverlapMeasureGenericFilter, LightObject);
+
+ //--------------------------------------------------------------------
+ void SetArgsInfo(const ArgsInfoType & a);
+ template<class FilterType>
+ void SetOptionsFromArgsInfoToFilter(FilterType * f) ;
+
+ //--------------------------------------------------------------------
+ // Main function called each time the filter is updated
+ template<class ImageType>
+ void UpdateWithInputImageType();
+
+ protected:
+ template<unsigned int Dim> void InitializeImageType();
+ ArgsInfoType mArgsInfo;
+
+ private:
+ LabelImageOverlapMeasureGenericFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+ };// end class
+ //--------------------------------------------------------------------
+} // end namespace clitk
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "clitkLabelImageOverlapMeasureGenericFilter.txx"
+#endif
+
+#endif // #define CLITKRELATIVEPOSITIONANALYZERGENERICFILTER_H
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+ ===========================================================================**/
+
+//--------------------------------------------------------------------
+template<class ArgsInfoType>
+clitk::LabelImageOverlapMeasureGenericFilter<ArgsInfoType>::
+LabelImageOverlapMeasureGenericFilter():
+ ImageToImageGenericFilter<Self>("LabelImageOverlapMeasure")
+{
+ // Default values
+ cmdline_parser_clitkDice_init(&mArgsInfo);
+ //InitializeImageType<2>();
+ InitializeImageType<3>();
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class ArgsInfoType>
+template<unsigned int Dim>
+void clitk::LabelImageOverlapMeasureGenericFilter<ArgsInfoType>::
+InitializeImageType()
+{
+ ADD_IMAGE_TYPE(Dim, uchar);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template<class ArgsInfoType>
+void clitk::LabelImageOverlapMeasureGenericFilter<ArgsInfoType>::
+SetArgsInfo(const ArgsInfoType & a)
+{
+ mArgsInfo=a;
+ //this->SetIOVerbose(mArgsInfo.verbose_flag);
+ if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
+ if (mArgsInfo.input1_given) this->AddInputFilename(mArgsInfo.input1_arg);
+ if (mArgsInfo.input2_given) this->AddInputFilename(mArgsInfo.input2_arg);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+// Update with the number of dimensions and the pixeltype
+//--------------------------------------------------------------------
+template<class ArgsInfoType>
+template<class FilterType>
+void clitk::LabelImageOverlapMeasureGenericFilter<ArgsInfoType>::
+SetOptionsFromArgsInfoToFilter(FilterType * f)
+{
+ f->SetLabel1(mArgsInfo.label1_arg);
+ f->SetLabel2(mArgsInfo.label2_arg);
+ f->SetVerboseFlag(mArgsInfo.verbose_flag);
+ f->SetLongFlag(mArgsInfo.long_flag);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+// Update with the number of dimensions and the pixeltype
+//--------------------------------------------------------------------
+template<class ArgsInfoType>
+template<class ImageType>
+void clitk::LabelImageOverlapMeasureGenericFilter<ArgsInfoType>::
+UpdateWithInputImageType()
+{
+ // Reading input
+ typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
+ typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
+
+ // Create filter
+ typedef clitk::LabelImageOverlapMeasureFilter<ImageType> FilterType;
+ typename FilterType::Pointer filter = FilterType::New();
+
+ // Set global Options
+ filter->SetInput(0, input1);
+ filter->SetInput(1, input2);
+ SetOptionsFromArgsInfoToFilter<FilterType>(filter);
+
+ // Go !
+ filter->Update();
+
+ // Write/Save results
+ // typename ImageType::Pointer output = filter->GetOutput();
+ // this->template SetNextOutput<ImageType>(output);
+}
+//--------------------------------------------------------------------
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+
+// clitk
+#include "clitkMatrixInverse_ggo.h"
+#include "clitkTransformUtilities.h"
+#include "clitkIO.h"
+#include "clitkMatrix.h"
+
+//--------------------------------------------------------------------
+int main(int argc, char * argv[])
+{
+ // Init command line
+ GGO(clitkMatrixInverse, args_info);
+ CLITK_INIT;
+
+ // Read matrix
+ itk::Matrix<double, 4, 4> matrix;
+ try {
+ matrix = clitk::ReadMatrix3D(args_info.input_arg);
+ }
+ catch (itk::ExceptionObject & err) {
+ std::cerr << "Error reading " << args_info.input_arg << std::endl;
+ std::cerr << err.GetDescription() << std::endl;
+ exit(-1);
+ }
+
+ matrix = matrix.GetInverse();
+
+ // Print matrix
+ std::ofstream os;
+ clitk::openFileForWriting(os, args_info.output_arg);
+ os << clitk::Get4x4MatrixDoubleAsString(matrix, 16);
+ os.close();
+
+ return EXIT_SUCCESS;
+}// end main
+
+//--------------------------------------------------------------------
--- /dev/null
+#File clitkMatrixToElastixTransform.ggo
+package "clitkMatrixInverse"
+version "1.0"
+purpose ""
+
+option "config" - "Config file" string optional
+option "input" i "Input matrix filename" string required
+option "output" o "Output matrix filename" string required
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+
+ Authors belong to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://www.centreleonberard.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+
+ - BSD See included LICENSE.txt file
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+===========================================================================**/
+
+// clitk
+#include "clitkMatrixToElastixTransform_ggo.h"
+#include "clitkTransformUtilities.h"
+#include "clitkIO.h"
+
+// itk
+#include <itkEuler3DTransform.h>
+
+//--------------------------------------------------------------------
+int main(int argc, char * argv[])
+{
+ // Init command line
+ GGO(clitkMatrixToElastixTransform, args_info);
+ CLITK_INIT;
+
+ // Read matrix
+ itk::Matrix<double, 4, 4> matrix;
+ try {
+ matrix = clitk::ReadMatrix3D(args_info.input_arg);
+ }
+ catch (itk::ExceptionObject & err) {
+ std::cerr << "Error reading " << args_info.input_arg << std::endl;
+ std::cerr << err.GetDescription() << std::endl;
+ exit(-1);
+ }
+
+ // Compute parameters from transfer using itk Euler transform
+ itk::Euler3DTransform<double>::CenterType center;
+ center.Fill(0.);
+ if(args_info.center_given==3) {
+ center[0] = args_info.center_arg[0];
+ center[1] = args_info.center_arg[1];
+ center[2] = args_info.center_arg[2];
+ }
+ itk::Euler3DTransform<double>::MatrixType rotMat;
+ itk::Euler3DTransform<double>::OutputVectorType transVec;
+ for(int i=0; i<3; i++) {
+ transVec[i] = matrix[i][3];
+ for(int j=0; j<3; j++)
+ rotMat[i][j] = matrix[i][j];
+ }
+ itk::Euler3DTransform<double>::Pointer euler;
+ euler = itk::Euler3DTransform<double>::New();
+ euler->SetCenter(center);
+ euler->SetOffset(transVec);
+ euler->SetComputeZYX(false);
+ try {
+ euler->SetMatrix(rotMat);
+ }
+ catch (itk::ExceptionObject & err) {
+ std::cerr << "Error reading " << args_info.input_arg << std::endl;
+ std::cerr << err.GetDescription() << std::endl;
+ exit(-1);
+ }
+
+ // Write result
+ std::ofstream out;
+ clitk::openFileForWriting(out, args_info.output_arg);
+ out << "(Transform \"EulerTransform\")" << std::endl;
+ out << "(NumberOfParameters 6)" << std::endl;
+ out << "(TransformParameters ";
+ for(unsigned int i=0; i<6; i++)
+ out << euler->GetParameters()[i] << ' ';
+ out << ')' << std::endl;
+ out << "(InitialTransformParametersFileName \"NoInitialTransform\")" << std::endl;
+ out << "(HowToCombineTransforms \"Compose\")" << std::endl;
+
+ out << "// EulerTransform specific" << std::endl;
+ out << "(CenterOfRotationPoint "<< center[0] << ' ' << center[1] << ' ' << center[2] << ')' << std::endl;
+ out << "(ComputeZYX \"false\")" << std::endl;
+
+ // The rest is commented, up to the user to define it manually
+ out << "// Image specific" << std::endl;
+ out << "// (FixedImageDimension 3)" << std::endl;
+ out << "// (MovingImageDimension 3)" << std::endl;
+ out << "// (FixedInternalImagePixelType \"float\")" << std::endl;
+ out << "// (MovingInternalImagePixelType \"float\")" << std::endl;
+ out << "// (Size 1 1 1)" << std::endl;
+ out << "// (Index 0 0 0)" << std::endl;
+ out << "// (Spacing 1 1 1)" << std::endl;
+ out << "// (Origin -0.5 -0.5 -0.5)" << std::endl;
+ out << "// (Direction 1 0 0 0 1 0 0 0 1)" << std::endl;
+ out << "// (UseDirectionCosines \"true\")" << std::endl << std::endl;
+
+ out << "// ResampleInterpolator specific" << std::endl;
+ out << "// (ResampleInterpolator \"FinalBSplineInterpolator\")" << std::endl;
+ out << "// (FinalBSplineInterpolationOrder 3)" << std::endl << std::endl;
+
+ out << "// Resampler specific" << std::endl;
+ out << "// (Resampler \"DefaultResampler\")" << std::endl;
+ out << "// (DefaultPixelValue 0.000000)" << std::endl;
+ out << "// (ResultImageFormat \"mhd\")" << std::endl;
+ out << "// (ResultImagePixelType \"short\")" << std::endl;
+ out << "// (CompressResultImage \"false\")" << std::endl;
+ out.close();
+
+ return EXIT_SUCCESS;
+}// end main
+
+//--------------------------------------------------------------------
--- /dev/null
+#File clitkMatrixToElastixTransform.ggo
+package "clitkMatrixToElastixTransform"
+version "1.0"
+purpose ""
+
+option "config" - "Config file" string optional
+option "verbose" v "Verbose" flag off
+option "input" i "Input matrix filename" string required
+option "output" o "Output elastix filename" string required
+option "center" - "Rotation center" float no multiple(3)
--- /dev/null
+/*=========================================================================
+ Program: vv http://www.creatis.insa-lyon.fr/rio/vv
+ Main authors : XX XX XX
+
+ Authors belongs to:
+ - University of LYON http://www.universite-lyon.fr/
+ - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
+ - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the copyright notices for more information.
+
+ It is distributed under dual licence
+ - BSD http://www.opensource.org/licenses/bsd-license.php
+ - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
+
+ =========================================================================*/
+
+#include <fstream>
+#include "clitkCommon.h"
+#include "clitkMergeAsciiDoseActor_ggo.h"
+
+//--------------------------------------------------------------------
+int main(int argc, char * argv[]) {
+
+ // Init command line
+ GGO(clitkMergeAsciiDoseActor, args_info);
+
+ // Read input 1
+ std::vector<double> input1;
+ clitk::readDoubleFromFile(args_info.input1_arg, input1);
+ DD(input1.size());
+
+ // Read input 2
+ std::vector<double> input2;
+ clitk::readDoubleFromFile(args_info.input2_arg, input2);
+ DD(input2.size());
+
+ // Check
+ if (input1.size() != input2.size()) {
+ std::cerr << "Error. Please provide input file with the same number of values. Read "
+ << input1.size() << " in " << args_info.input1_arg
+ << " and " << input2.size() << " in " << args_info.input2_arg << std::endl;
+ exit(0);
+ }
+
+ // Add
+ for(unsigned int i=0; i<input1.size(); i++) {
+ input1[i] += input2[i];
+ }
+
+ // Write output
+ std::ofstream os;
+ clitk::openFileForWriting(os, args_info.output_arg);
+ os << "## Merge " << input1.size() << " values" << std::endl;
+ for(unsigned int i=0; i<input1.size(); i++) {
+ os << input1[i] << std::endl;
+ }
+ os.close();
+
+ // This is the end my friend
+ return 0;
+}
--- /dev/null
+# file clitkMergeAsciiDoseActor.ggo
+package "clitk"
+version "Add values from two ASCII files (must contains list of numbers)"
+
+option "config" - "Config file" string no
+option "input1" i "Input ASCII file" string yes
+option "input2" j "Input ASCII file" string yes
+option "output" o "Output ASCII file" string yes
+
+
--- /dev/null
+/**
+ =================================================
+ * @file clitkMergeRootFiles.cxx
+ * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
+ * @date 01 Apr 2009
+ *
+ * @brief
+ *
+ =================================================*/
+
+#include "clitkMergeRootFiles_ggo.h"
+#include "clitkCommon.h"
+#include "GateMergeManager.hh"
+#include <string>
+#include <TROOT.h>
+#include <TPluginManager.h>
+#include <TFile.h>
+#include <TKey.h>
+#include <TFileMerger.h>
+#include <TTree.h>
+#include <TChain.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <iostream>
+
+using std::endl;
+using std::cout;
+
+struct PetInputFile
+{
+ string filename;
+ double mean_time;
+};
+
+bool sort_pet_input_file(const PetInputFile& a, const PetInputFile& b)
+{
+ return a.mean_time<b.mean_time;
+};
+
+typedef std::vector<PetInputFile> PetInputFiles;
+
+//-----------------------------------------------------------------------------
+int main(int argc, char * argv[]) {
+
+ gROOT->GetPluginManager()->AddHandler("TVirtualStreamerInfo", "*",
+ "TStreamerInfo", "RIO", "TStreamerInfo()");
+
+ // init command line
+ GGO(clitkMergeRootFiles, args_info);
+
+ // Check parameters
+ if (args_info.input_given < 2) {
+ FATAL("Error, please provide at least two inputs files");
+ }
+
+ { // Detect Pet output
+ bool all_pet_output = true;
+ PetInputFiles pet_input_files;
+ for (uint i=0; i<args_info.input_given; i++)
+ {
+ const char* filename = args_info.input_arg[i];
+ PetInputFile input_file;
+ input_file.filename = filename;
+ TFile* handle = TFile::Open(filename,"READ");
+ assert(handle);
+ TTree* hits = dynamic_cast<TTree*>(handle->Get("Hits"));
+ TTree* singles = dynamic_cast<TTree*>(handle->Get("Singles"));
+ const bool is_pet_output = (hits!=NULL) && (singles!=NULL);
+ cout << "testing " << filename << " is_pet_output " << is_pet_output;
+
+ if (is_pet_output)
+ {
+ double time;
+ double time_accum = 0;
+ singles->SetBranchAddress("time",&time);
+ size_t total = singles->GetEntries();
+ for (size_t kk=0; kk<total; kk++)
+ {
+ singles->GetEntry(kk);
+ time_accum += time;
+ }
+
+ input_file.mean_time = time_accum/total;
+ pet_input_files.push_back(input_file);
+ cout << " mean_time " << input_file.mean_time;
+ }
+
+ cout << endl;
+
+ handle->Close();
+ delete handle;
+ all_pet_output &= is_pet_output;
+ }
+ cout << "all_pet_output " << all_pet_output << endl;
+
+ if (all_pet_output)
+ {
+ GateMergeManager manager(args_info.fastmerge_given,args_info.verbose_arg,true,0,"");
+
+ cout << "sorting input file using singles time" << endl;
+ std::sort(pet_input_files.begin(),pet_input_files.end(),sort_pet_input_file);
+
+ Strings input_filenames;
+ for (PetInputFiles::const_iterator iter=pet_input_files.begin(); iter!=pet_input_files.end(); iter++)
+ input_filenames.push_back(iter->filename);
+
+ manager.StartMergingFromFilenames(input_filenames,args_info.output_arg);
+ return 0;
+ }
+ }
+
+
+ // Merge
+ TFileMerger * merger = new TFileMerger;
+ for (uint i=0; i<args_info.input_given; i++) merger->AddFile(args_info.input_arg[i]);
+ merger->OutputFile(args_info.output_arg);
+ merger->Merge();
+
+ // this is the end my friend
+ return 0;
+}
+//-----------------------------------------------------------------------------
--- /dev/null
+# --------------------------------------------
+# file clitkMergeRootFiles.ggo
+package "clitkMergeRootFiles"
+version "1.0"
+purpose "Merge several ROOT files"
+# --------------------------------------------
+
+option "config" - "Config file" string no
+option "input" i "Input ROOT filenames" string multiple yes
+option "output" o "Output ROOT filename" string yes
+option "fastmerge" f "Fast merge" optional
+option "verbose" v "Verbose level" int optional
+
#include "vvSaveState.h"
#include "vvReadState.h"
#include "clitkConfiguration.h"
+#include "clitkMatrix.h"
// ITK include
#include <itkImage.h>
infoPanel->setNPixel(QString::number(NPixel)+" ("+inputSizeInBytes+")");
transformation = imageSelected->GetTransform()[tSlice]->GetMatrix();
- infoPanel->setTransformation(Get4x4MatrixDoubleAsString(transformation));
+ infoPanel->setTransformation(clitk::Get4x4MatrixDoubleAsString(transformation).c_str());
landmarksPanel->SetCurrentLandmarks(mSlicerManagers[index]->GetLandmarks(),
mSlicerManagers[index]->GetTSlice());
}
//------------------------------------------------------------------------------
-//------------------------------------------------------------------------------
-QString vvMainWindow::Get4x4MatrixDoubleAsString(vtkSmartPointer<vtkMatrix4x4> matrix, const int precision)
-{
- std::ostringstream strmatrix;
-
- // Figure out the number of digits of the integer part of the largest absolute value
- // for each column
- unsigned width[4];
- for (unsigned int j = 0; j < 4; j++){
- double absmax = 0.;
- for (unsigned int i = 0; i < 4; i++)
- absmax = std::max(absmax, vnl_math_abs(matrix->GetElement(i, j)));
- unsigned ndigits = (unsigned)std::max(0.,std::log10(absmax))+1;
- width[j] = precision+ndigits+3;
- }
-
- // Output with correct width, aligned to the right
- for (unsigned int i = 0; i < 4; i++) {
- for (unsigned int j = 0; j < 4; j++) {
- strmatrix.setf(ios::fixed,ios::floatfield);
- strmatrix.precision(precision);
- strmatrix.fill(' ');
- strmatrix.width(width[j]);
- strmatrix << std::right << matrix->GetElement(i, j);
- }
- strmatrix << std::endl;
- }
- QString result = strmatrix.str().c_str();
- return result;
-}
-//------------------------------------------------------------------------------
-
//------------------------------------------------------------------------------
QString vvMainWindow::GetVectorDoubleAsString(std::vector<double> vectorDouble)
{
void SaveCurrentStateAs(const std::string& stateFile);
void ReadSavedStateFile(const std::string& stateFile);
void LinkAllImages();
- QString Get4x4MatrixDoubleAsString(vtkSmartPointer<vtkMatrix4x4> matrix, const int precision=3);
virtual void UpdateCurrentSlicer();
virtual QTabWidget * GetTab();
// clitk
#include "clitkTransformUtilities.h"
+#include "clitkMatrix.h"
// qt
#include <QMessageBox>
for(int i=0; i<4; i++)
// TODO SR and BP: check on the list of transforms and not the first only
mInitialMatrix->SetElement(i,j, mCurrentSlicerManager->GetImage()->GetTransform()[0]->GetMatrix()->GetElement(i,j));
- QString origTransformString = dynamic_cast<vvMainWindow*>(mMainWindow)->Get4x4MatrixDoubleAsString(mInitialMatrix);
+ QString origTransformString(clitk::Get4x4MatrixDoubleAsString(mInitialMatrix).c_str());
transformationLabel->setText(origTransformString);
SetTransform(mInitialMatrix);
if (file.open(QFile::WriteOnly | QFile::Truncate)) {
// TODO SR and BP: check on the list of transforms and not the first only
vtkMatrix4x4* matrix = mCurrentSlicerManager->GetImage()->GetTransform()[0]->GetMatrix();
- QString matrixStr = dynamic_cast<vvMainWindow*>(mMainWindow)->Get4x4MatrixDoubleAsString(matrix,16);
+ QString matrixStr = clitk::Get4x4MatrixDoubleAsString(matrix,16).c_str();
QTextStream out(&file);
out << matrixStr;
}
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff