Skip to content

Commit

Permalink
BUG: Spatial Object Registration Optimization (#86)
Browse files Browse the repository at this point in the history 
BUG: Ensure consistent metric values for SpatialObject registration optimizations
BUG: Remove github workflow for Python3.6.
BUG: Also update linux wheel name to _x86_64
  • Loading branch information
@aylward
aylward authored Oct 7, 2021
1 parent 119f23c commit 3911cc7
Show file tree
Hide file tree
Showing 7 changed files with 146 additions and 131 deletions.
12 changes: 7 additions & 5 deletions .github/workflows/build-test-package.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ jobs:
sudo apt-get install -y libgl1-mesa-glx
sudo apt-get install -y libglvnd-core-dev
sudo apt-get install -y freeglut3-dev
sudo apt-get install -y libtbb-dev
- name: Install build dependencies
run: |
Expand Down Expand Up @@ -201,7 +202,7 @@ jobs:
strategy:
max-parallel: 1
matrix:
python-version: [36, 37, 38, 39]
python-version: [37, 38, 39]
include:
- itk-python-git-tag: "v5.3rc01"
c-compiler: "gcc"
Expand Down Expand Up @@ -229,6 +230,7 @@ jobs:
sudo apt-get install -y libgl1-mesa-glx
sudo apt-get install -y libglvnd-core-dev
sudo apt-get install -y freeglut3-dev
sudo apt-get install -y libtbb-dev
python -m pip install --upgrade pip
python -m pip install ninja
Expand All @@ -254,16 +256,16 @@ jobs:
run: |
cat > git-configure-build-sem.sh << EOF
#!/bin/sh
ln -s -f -T /ITKPythonPackage/ITK-cp36-cp36m-manylinux2014_x64 /work/ITK-cp36-cp36m-manylinux2014_x64
ln -s -f -T /ITKPythonPackage/ITK-cp37-cp37m-manylinux2014_x86_64 /work/ITK-cp37-cp37m-manylinux2014_x86_64
ln -s -f -T /ITKPythonPackage/ITK-source /work/ITK-source
export ITK_DIR="/work/ITK-cp36-cp36m-manylinux2014_x64"
export ITK_DIR="/work/ITK-cp37-cp37m-manylinux2014_x86_64"
echo "itk_dir = ${ITK_DIR}"
ls -la /work
git clone https://github.com/Slicer/SlicerExecutionModel
mkdir -p SlicerExecutionModel-build
cd SlicerExecutionModel-build
cmake -DCMAKE_BUILD_TYPE:STRING=Release \
-DITK_DIR=/work/ITK-cp36-cp36m-manylinux2014_x64 \
-DITK_DIR=/work/ITK-cp37-cp37m-manylinux2014_x86_64 \
-GNinja \
../SlicerExecutionModel
ninja
Expand Down Expand Up @@ -393,7 +395,7 @@ jobs:
strategy:
max-parallel: 1
matrix:
python-version-minor: [6, 7, 8, 9]
python-version-minor: [7, 8, 9]
include:
- itk-python-git-tag: "v5.3rc01"
c-compiler: "cl.exe"
Expand Down
10 changes: 4 additions & 6 deletions src/Numerics/itktubeNJetImageFunction.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -592,7 +592,7 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,
// EVALUATE
if( m_UseProjection )
{
return EvaluateAtContinuousIndex( cIndex, scale );
m_MostRecentIntensity = EvaluateAtContinuousIndex( cIndex, scale );
}
else
{
Expand All @@ -615,9 +615,8 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,

m_MostRecentIntensity = ( val0 + 0.5455*val1 + 0.5455*val2 )
/ ( 1+2*0.5455 ) / 2;

return m_MostRecentIntensity;
}
return m_MostRecentIntensity;
}


Expand All @@ -630,7 +629,7 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,
// EVALUATE
if( m_UseProjection )
{
return EvaluateAtContinuousIndex( cIndex, scale );
m_MostRecentIntensity = EvaluateAtContinuousIndex( cIndex, scale );
}
else
{
Expand Down Expand Up @@ -663,9 +662,8 @@ EvaluateAtContinuousIndex( const ContinuousIndexType & cIndex,

m_MostRecentIntensity =
( 2*val0 + 0.5455*( val1+val2+val3+val4 ) )/( 2+4*0.5455 ) / 2;

return m_MostRecentIntensity;
}
return m_MostRecentIntensity;
}

template< class TInputImage >
Expand Down
204 changes: 113 additions & 91 deletions src/Registration/itktubeOptimizedSpatialObjectToImageRegistrationMethod.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,6 @@ void
OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
::GenerateData( void )
{
std::cout << "ORM: GenerateData" << std::endl;
if( this->GetReportProgress() )
{
std::cout << "UPDATE START" << std::endl;
Expand Down Expand Up @@ -225,10 +224,8 @@ OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
this->GetMovingSpatialObjectMaskObject() );
}

std::cout << "ORM: GenerateData: Metric Init" << std::endl;
metric->Initialize();

std::cout << "ORM: GenerateData: Optimize" << std::endl;
try
{
this->Optimize(metric);
Expand All @@ -250,100 +247,113 @@ void
OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
::Optimize( MetricType * metric )
{
if( m_UseEvolutionaryOptimization && this->GetMaxIterations() > 0 )
if( this->GetMaxIterations() > 0 )
{
if( this->GetReportProgress() )
{
std::cout << "EVOLUTIONARY START" << std::endl;
}

typedef OnePlusOneEvolutionaryOptimizer EvoOptimizerType;
EvoOptimizerType::Pointer evoOpt = EvoOptimizerType::New();

evoOpt->SetNormalVariateGenerator( Statistics::NormalVariateGenerator
::New() );
evoOpt->SetEpsilon( this->GetTargetError() );
evoOpt->Initialize( 0.1 );
evoOpt->SetCatchGetValueException( true );
evoOpt->SetMetricWorstPossibleValue( 9999999 );
evoOpt->SetScales( this->GetTransformParametersScales() );
evoOpt->SetMaximumIteration( this->GetMaxIterations() );
evoOpt->SetMaximize( false );

if( this->GetObserver() )
{
evoOpt->AddObserver( IterationEvent(), this->GetObserver() );
}

if( this->GetReportProgress() )
if( m_UseEvolutionaryOptimization )
{
typedef SpatialObjectToImageRegistrationViewer ViewerCommandType;
typename ViewerCommandType::Pointer command = ViewerCommandType::New();
if( this->GetTransform()->GetNumberOfParameters() > 16 )
if( this->GetReportProgress() )
{
command->SetDontShowParameters( true );
std::cout << "EVOLUTIONARY START" << std::endl;
}

typedef OnePlusOneEvolutionaryOptimizer EvoOptimizerType;
EvoOptimizerType::Pointer evoOpt = EvoOptimizerType::New();

evoOpt->SetNormalVariateGenerator( Statistics::NormalVariateGenerator
::New() );
evoOpt->SetEpsilon( this->GetTargetError() );
evoOpt->Initialize( 0.1 );
evoOpt->SetCatchGetValueException( true );
evoOpt->SetMetricWorstPossibleValue( 0 );
evoOpt->SetScales( this->GetTransformParametersScales() );
evoOpt->SetMaximumIteration( this->GetMaxIterations() * 0.75 );

if( this->GetObserver() )
{
evoOpt->AddObserver( IterationEvent(), this->GetObserver() );
}

if( this->GetReportProgress() )
{
typedef SpatialObjectToImageRegistrationViewer ViewerCommandType;
typename ViewerCommandType::Pointer command = ViewerCommandType::New();
if( this->GetTransform()->GetNumberOfParameters() > 16 )
{
command->SetDontShowParameters( true );
}
evoOpt->AddObserver( IterationEvent(), command );
}

evoOpt->SetCostFunction(metric);
evoOpt->SetInitialPosition(m_InitialTransformParameters);
try
{
evoOpt->StartOptimization();
}
catch (...)
{
std::cout << "Exception caught during evolutionary registration."
<< std::endl;
std::cout << "Continuing using best values..." << std::endl;
std::cout << " Pos = " << evoOpt->GetCurrentPosition()
<< std::endl << std::endl;
}

if( evoOpt->GetCurrentPosition().size() ==
this->GetTransform()->GetNumberOfParameters() )
{
bool valid = true;
for(unsigned int i=0; i<evoOpt->GetCurrentPosition().size(); ++i)
{
if(evoOpt->GetCurrentPosition()[i] != evoOpt->GetCurrentPosition()[i])
{
valid = false;
break;
}
}
if( valid )
{
this->SetLastTransformParameters( evoOpt->GetCurrentPosition() );
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
}
else
{
std::cerr << "Error: Invalid Evolutionalry final parameters"
<< std::endl;
this->SetLastTransformParameters( m_InitialTransformParameters );
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = 0;
}
}
if( this->GetReportProgress() )
{
std::cout << "EVOLUTIONARY END" << std::endl;
}
evoOpt->AddObserver( IterationEvent(), command );
}

evoOpt->SetCostFunction(metric);
evoOpt->SetInitialPosition(m_InitialTransformParameters);
try
{
evoOpt->StartOptimization();
}
catch (...)
{
std::cout << "Exception caught during evolutionary registration."
<< std::endl;
std::cout << "Continuing using best values..." << std::endl;
std::cout << " Pos = " << evoOpt->GetCurrentPosition()
<< std::endl << std::endl;
}

if( evoOpt->GetCurrentPosition().size() !=
this->GetTransform()->GetNumberOfParameters() )
{
this->SetLastTransformParameters( this->GetInitialTransformParameters() );
}
else
{
this->SetLastTransformParameters( evoOpt->GetCurrentPosition() );
this->SetLastTransformParameters( m_InitialTransformParameters );
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = 0;
}
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
}
else
{
this->SetLastTransformParameters( this->GetInitialTransformParameters() );
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = 0;
}

if( this->GetReportProgress() )
{
std::cout << "EVOLUTIONARY END" << std::endl;
}

if( this->GetMaxIterations() > 0 )
{
//if( this->GetReportProgress() )
if( this->GetReportProgress() )
{
std::cout << "GRADIENT START" << std::endl;
}

typedef FRPROptimizer GradOptimizerType;
GradOptimizerType::Pointer gradOpt = GradOptimizerType::New();

gradOpt->SetMaximize( false );
gradOpt->SetCatchGetValueException( true );
gradOpt->SetMetricWorstPossibleValue( 9999999 );
gradOpt->SetStepLength( 0.25 );
gradOpt->SetMetricWorstPossibleValue( 0 );
gradOpt->SetStepLength( 0.1 );
gradOpt->SetStepTolerance( this->GetTargetError() );
gradOpt->SetMaximumIteration( this->GetMaxIterations() );
gradOpt->SetMaximumIteration( this->GetMaxIterations() * 0.25 / 10 );
gradOpt->SetMaximumLineIteration( 10 );
gradOpt->SetScales( this->GetTransformParametersScales() );
gradOpt->SetUseUnitLengthGradient(true);
gradOpt->SetUseUnitLengthGradient(false);
gradOpt->SetToFletchReeves();

//if( this->GetReportProgress() )
Expand All @@ -362,7 +372,7 @@ OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
}

gradOpt->SetCostFunction(metric);
gradOpt->SetInitialPosition(m_InitialTransformParameters);
gradOpt->SetInitialPosition(m_LastTransformParameters);
try
{
gradOpt->StartOptimization();
Expand All @@ -376,30 +386,42 @@ OptimizedSpatialObjectToImageRegistrationMethod<ObjectDimension, TImage>
<< std::endl << std::endl;
}

if( gradOpt->GetCurrentPosition().size() !=
if( gradOpt->GetCurrentPosition().size() ==
this->GetTransform()->GetNumberOfParameters() )
{
this->SetLastTransformParameters( this->GetInitialTransformParameters() );
bool valid = true;
for(unsigned int i=0; i<gradOpt->GetCurrentPosition().size(); ++i)
{
if(gradOpt->GetCurrentPosition()[i] != gradOpt->GetCurrentPosition()[i])
{
valid = false;
break;
}
}
if( valid )
{
this->SetLastTransformParameters( gradOpt->GetCurrentPosition() );
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
}
else
{
this->SetLastTransformParameters( m_InitialTransformParameters );
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = 0;
}
}
else
if( this->GetReportProgress() )
{
this->SetLastTransformParameters( gradOpt->GetCurrentPosition() );
std::cout << "GRADIENT END" << std::endl;
}
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = metric->GetValue( m_LastTransformParameters );
}
else
{
this->SetLastTransformParameters( this->GetInitialTransformParameters() );
this->SetLastTransformParameters( m_InitialTransformParameters );
this->GetTransform()->SetParametersByValue( m_LastTransformParameters );
m_FinalMetricValue = 0;
}


if( this->GetReportProgress() )
{
std::cout << "GRADIENT END" << std::endl;
}
}

template <unsigned int ObjectDimension, class TImage>
Expand Down
Loading

0 comments on commit 3911cc7

Please sign in to comment.