<commit_before>/*=========================================================================

Library:   TubeTK

Copyright 2010 Kitware Inc. 28 Corporate Drive,
Clifton Park, NY, 12065, USA.

All rights reserved.

Licensed under the Apache License, Version 2.0 ( the "License" );
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

=========================================================================*/

#include "itktubeRidgeExtractor.h"

#include <itkImageFileReader.h>
#include <itkMersenneTwisterRandomVariateGenerator.h>
#include <itkSpatialObjectReader.h>

int itktubeRidgeExtractorTest2( int argc, char * argv[] )
{
  if( argc != 3 )
    {
    std::cout
      << "itktubeRidgeExtractorTest <inputImage> <vessel.tre>"
      << std::endl;
    return EXIT_FAILURE;
    }

  typedef itk::Image<float, 3>   ImageType;

  typedef itk::ImageFileReader< ImageType > ImageReaderType;
  ImageReaderType::Pointer imReader = ImageReaderType::New();
  imReader->SetFileName( argv[1] );
  imReader->Update();

  ImageType::Pointer im = imReader->GetOutput();

  typedef itk::tube::RidgeExtractor<ImageType> RidgeOpType;
  RidgeOpType::Pointer ridgeOp = RidgeOpType::New();

  //ridgeOp->SetDebug( true );

  ridgeOp->SetInputImage( im );
  ridgeOp->SetStepX( 0.5 );
  ridgeOp->SetDynamicScale( true );

  typedef itk::SpatialObjectReader<>                   ReaderType;
  typedef itk::SpatialObject<>::ChildrenListType       ObjectListType;
  typedef itk::GroupSpatialObject<>                    GroupType;
  typedef itk::TubeSpatialObject<>                     TubeType;
  typedef TubeType::TubePointListType                  TubePointListType;
  typedef TubeType::TubePointType                      TubePointType;

  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( argv[2] );
  reader->Update();
  GroupType::Pointer group = reader->GetGroup();

  std::cout << "Number of children = " << group->GetNumberOfChildren()
    << std::endl;

  char tubeName[17];
  std::strcpy( tubeName, "Tube" );
  ObjectListType * tubeList = group->GetChildren( -1, tubeName );

  unsigned int numTubes = tubeList->size();
  std::cout << "Number of tubes = " << numTubes << std::endl;

  typedef itk::Statistics::MersenneTwisterRandomVariateGenerator
    RandGenType;
  RandGenType::Pointer rndGen = RandGenType::New();
  rndGen->Initialize(); // set seed here

  RidgeOpType::IndexType imMinI = ridgeOp->GetExtractBoundMinInIndexSpace();
  RidgeOpType::IndexType imMaxI = ridgeOp->GetExtractBoundMaxInIndexSpace();
  int margin = 3;
  std::cout << "Bound min = " << imMinI << std::endl;
  std::cout << "Bound max = " << imMaxI << std::endl;

  int failures = 0;
  for( int mcRun=0; mcRun<2; mcRun++ )
    {
    std::cout << std::endl;
    std::cout << std::endl;
    std::cout << "***** Beginning tube test ***** " << std::endl;
    unsigned int rndTubeNum = rndGen->GetUniformVariate( 0, 1 ) * numTubes;
    if( rndTubeNum > numTubes-1 )
      {
      rndTubeNum = numTubes-1;
      }
    ObjectListType::iterator tubeIter = tubeList->begin();
    for( unsigned int i=0; i<rndTubeNum; i++ )
      {
      ++tubeIter;
      }
    TubeType::Pointer tube = static_cast< TubeType * >(
      tubeIter->GetPointer() );
    std::cout << "Test tube = " << rndTubeNum << std::endl;

    TubePointListType tubePointList = tube->GetPoints();
    unsigned int numPoints = tubePointList.size();
    unsigned int rndPointNum = rndGen->GetUniformVariate( 0, 1 )
      * numPoints;
    if( rndPointNum > numPoints-1 )
      {
      rndPointNum = numPoints-1;
      }
    TubePointListType::iterator pntIter = tubePointList.begin();
    for( unsigned int i=0; i<rndPointNum; i++ )
      {
      ++pntIter;
      }
    TubePointType * pnt = static_cast< TubePointType * >( &( *pntIter ) );
    std::cout << "Test point = " << rndPointNum << std::endl;

    tube->Update();

    ImageType::PointType pntX = pnt->GetPositionInWorldSpace();
    std::cout << "Test position = " << pntX << std::endl;

    RidgeOpType::ContinuousIndexType xContI;
    bool inMargin = im->TransformPhysicalPointToContinuousIndex( pntX,
      xContI );
    if( inMargin )
      {
      for( unsigned int i=0; i<ImageType::ImageDimension; i++ )
        {
        if( xContI[i] < imMinI[i]+margin || xContI[i] > imMaxI[i]-margin )
          {
          inMargin = false;
          break;
          }
        }
      }
    if( !inMargin )
      {
      --mcRun;
      std::cout << "Warning: Tube point outside of image, repicking..."
        << std::endl;
      std::cout << "   Point Index = " << xContI << std::endl;
      std::cout << "   Min Index = " << imMinI << std::endl;
      std::cout << "   Max Index = " << imMaxI << std::endl;
      std::cout << "   Margin = " << margin << std::endl;
      continue;
      }

    RidgeOpType::IndexType minI = imMinI;
    RidgeOpType::IndexType maxI = imMaxI;
    if( xContI[2] - margin > minI[2] )
      {
      minI[2] = xContI[2] - margin;
      }
    if( xContI[2] + margin < maxI[2] )
      {
      maxI[2] = xContI[2] + margin;
      }
    std::cout << "Setting min and max z to save time." << std::endl;
    std::cout << "   Min index = " << minI << std::endl;
    std::cout << "   Max index = " << maxI << std::endl;
    ridgeOp->SetExtractBoundMinInIndexSpace( minI );
    ridgeOp->SetExtractBoundMaxInIndexSpace( maxI );

    if( pnt->GetRadiusInObjectSpace() > 1 )
      {
      ridgeOp->SetScale( 0.8 * pnt->GetRadiusInObjectSpace() );
      }
    else
      {
      ridgeOp->SetScale( 0.5 );
      }
    std::cout << "Initial scale = " << ridgeOp->GetScale() << std::endl;


    RidgeOpType::PointType xRidgePnt;
    im->TransformContinuousIndexToPhysicalPoint( xContI, xRidgePnt );
    if( ridgeOp->LocalRidge( xRidgePnt ) != RidgeOpType::SUCCESS )
      {
      RidgeOpType::ContinuousIndexType xRidgeContI;
      im->TransformPhysicalPointToContinuousIndex( xRidgePnt, xRidgeContI );
      std::cout << "*** FAILURE: Local ridge test failed.  No ridge found."
        << std::endl;
      std::cout << "   Source = " << xContI << std::endl;
      std::cout << "   Result = " << xRidgeContI << std::endl;
      ++failures;
      continue;
      }
    RidgeOpType::ContinuousIndexType xRidgeContI;
    im->TransformPhysicalPointToContinuousIndex( xRidgePnt, xRidgeContI );

    double diff = 0;
    for( unsigned int i=0; i<ImageType::ImageDimension; i++ )
      {
      double tf = xContI[i]-xRidgeContI[i];
      diff += tf * tf;
      }
    diff = std::sqrt( diff );
    if( diff > 2*pnt->GetRadiusInObjectSpace() && diff > 4 )
      {
      std::cout << "*** FAILURE: Local ridge test failed.  Local ridge too far."
        << std::endl;
      std::cout << "   Source = " << xContI << std::endl;
      std::cout << "   Result = " << xRidgeContI << std::endl;
      ++failures;
      continue;
      }

    std::cout << "Local ridge discovery a success!" << std::endl;
    std::cout << std::endl;
    std::cout << "***** Beginning tube extraction ***** " << std::endl;<commit_msg>COMP: Add more debugging to test that intermittently fails<commit_after>/*=========================================================================

Library:   TubeTK

Copyright 2010 Kitware Inc. 28 Corporate Drive,
Clifton Park, NY, 12065, USA.

All rights reserved.

Licensed under the Apache License, Version 2.0 ( the "License" );
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

=========================================================================*/

#include "itktubeRidgeExtractor.h"

#include <itkImageFileReader.h>
#include <itkMersenneTwisterRandomVariateGenerator.h>
#include <itkSpatialObjectReader.h>

int itktubeRidgeExtractorTest2( int argc, char * argv[] )
{
  if( argc != 3 )
    {
    std::cout
      << "itktubeRidgeExtractorTest <inputImage> <vessel.tre>"
      << std::endl;
    return EXIT_FAILURE;
    }

  typedef itk::Image<float, 3>   ImageType;

  typedef itk::ImageFileReader< ImageType > ImageReaderType;
  ImageReaderType::Pointer imReader = ImageReaderType::New();
  imReader->SetFileName( argv[1] );
  imReader->Update();

  ImageType::Pointer im = imReader->GetOutput();

  typedef itk::tube::RidgeExtractor<ImageType> RidgeOpType;
  RidgeOpType::Pointer ridgeOp = RidgeOpType::New();

  //ridgeOp->SetDebug( true );

  ridgeOp->SetInputImage( im );
  ridgeOp->SetStepX( 0.5 );
  ridgeOp->SetDynamicScale( true );

  typedef itk::SpatialObjectReader<>                   ReaderType;
  typedef itk::SpatialObject<>::ChildrenListType       ObjectListType;
  typedef itk::GroupSpatialObject<>                    GroupType;
  typedef itk::TubeSpatialObject<>                     TubeType;
  typedef TubeType::TubePointListType                  TubePointListType;
  typedef TubeType::TubePointType                      TubePointType;

  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( argv[2] );
  reader->Update();
  GroupType::Pointer group = reader->GetGroup();

  std::cout << "Number of children = " << group->GetNumberOfChildren()
    << std::endl;

  char tubeName[17];
  std::strcpy( tubeName, "Tube" );
  ObjectListType * tubeList = group->GetChildren( -1, tubeName );

  unsigned int numTubes = tubeList->size();
  std::cout << "Number of tubes = " << numTubes << std::endl;

  typedef itk::Statistics::MersenneTwisterRandomVariateGenerator
    RandGenType;
  RandGenType::Pointer rndGen = RandGenType::New();
  rndGen->Initialize(); // set seed here

  RidgeOpType::IndexType imMinI = ridgeOp->GetExtractBoundMinInIndexSpace();
  RidgeOpType::IndexType imMaxI = ridgeOp->GetExtractBoundMaxInIndexSpace();
  int margin = 3;
  std::cout << "Bound min = " << imMinI << std::endl;
  std::cout << "Bound max = " << imMaxI << std::endl;

  int failures = 0;
  for( int mcRun=0; mcRun<2; mcRun++ )
    {
    std::cout << std::endl;
    std::cout << std::endl;
    std::cout << "***** Beginning tube test ***** " << std::endl;
    unsigned int rndTubeNum = rndGen->GetUniformVariate( 0, 1 ) * numTubes;
    if( rndTubeNum > numTubes-1 )
      {
      rndTubeNum = numTubes-1;
      }
    ObjectListType::iterator tubeIter = tubeList->begin();
    for( unsigned int i=0; i<rndTubeNum; i++ )
      {
      ++tubeIter;
      }
    TubeType::Pointer tube = static_cast< TubeType * >(
      tubeIter->GetPointer() );
    std::cout << "Test tube = " << rndTubeNum << std::endl;

    TubePointListType tubePointList = tube->GetPoints();
    unsigned int numPoints = tubePointList.size();
    unsigned int rndPointNum = rndGen->GetUniformVariate( 0, 1 )
      * numPoints;
    if( rndPointNum > numPoints-1 )
      {
      rndPointNum = numPoints-1;
      }
    TubePointListType::iterator pntIter = tubePointList.begin();
    for( unsigned int i=0; i<rndPointNum; i++ )
      {
      ++pntIter;
      }
    TubePointType * pnt = static_cast< TubePointType * >( &( *pntIter ) );
    std::cout << "Test point = " << rndPointNum << std::endl;

    tube->Update();

    ImageType::PointType pntX = pnt->GetPositionInWorldSpace();
    std::cout << "Test position = " << pntX << std::endl;

    RidgeOpType::ContinuousIndexType xContI;
    bool inMargin = im->TransformPhysicalPointToContinuousIndex( pntX,
      xContI );
    if( inMargin )
      {
      for( unsigned int i=0; i<ImageType::ImageDimension; i++ )
        {
        if( xContI[i] < imMinI[i]+margin || xContI[i] > imMaxI[i]-margin )
          {
          inMargin = false;
          break;
          }
        }
      }
    if( !inMargin )
      {
      --mcRun;
      std::cout << "Warning: Tube point outside of image, repicking..."
        << std::endl;
      std::cout << "   Point Index = " << xContI << std::endl;
      std::cout << "   Min Index = " << imMinI << std::endl;
      std::cout << "   Max Index = " << imMaxI << std::endl;
      std::cout << "   Margin = " << margin << std::endl;
      continue;
      }

    RidgeOpType::IndexType minI = imMinI;
    RidgeOpType::IndexType maxI = imMaxI;
    if( xContI[2] - margin > minI[2] )
      {
      minI[2] = xContI[2] - margin;
      }
    if( xContI[2] + margin < maxI[2] )
      {
      maxI[2] = xContI[2] + margin;
      }
    std::cout << "Setting min and max z to save time." << std::endl;
    std::cout << "   Min index = " << minI << std::endl;
    std::cout << "   Max index = " << maxI << std::endl;
    ridgeOp->SetExtractBoundMinInIndexSpace( minI );
    ridgeOp->SetExtractBoundMaxInIndexSpace( maxI );

    if( pnt->GetRadiusInObjectSpace() > 1 )
      {
      ridgeOp->SetScale( 0.8 * pnt->GetRadiusInObjectSpace() );
      }
    else
      {
      ridgeOp->SetScale( 0.5 );
      }
    std::cout << "Initial scale = " << ridgeOp->GetScale() << std::endl;


    RidgeOpType::PointType xRidgePnt;
    im->TransformContinuousIndexToPhysicalPoint( xContI, xRidgePnt );
    if( ridgeOp->LocalRidge( xRidgePnt ) != RidgeOpType::SUCCESS )
      {
      RidgeOpType::ContinuousIndexType xRidgeContI;
      im->TransformPhysicalPointToContinuousIndex( xRidgePnt, xRidgeContI );
      std::cout << "*** FAILURE: Local ridge test failed.  No ridge found."
        << std::endl;
      std::cout << "   Source = " << xContI << std::endl;
      std::cout << "   Result = " << xRidgeContI << std::endl;
      double ridgeness, intensity, roundness, curvature, levelness;
      ridgeness = ridgeOp->Ridgeness( xRidgePnt, intensity, roundness,
        curvature, levelness );
      std::cout << "       ridgness = " << ridgeness << std::endl;
      std::cout << "       intensity = " << intensity << std::endl;
      std::cout << "       roundness = " << roundness << std::endl;
      std::cout << "       curvature = " << curvature << std::endl;
      std::cout << "       levelness = " << levelness << std::endl;
      ++failures;
      continue;
      }
    RidgeOpType::ContinuousIndexType xRidgeContI;
    im->TransformPhysicalPointToContinuousIndex( xRidgePnt, xRidgeContI );

    double diff = 0;
    for( unsigned int i=0; i<ImageType::ImageDimension; i++ )
      {
      double tf = xContI[i]-xRidgeContI[i];
      diff += tf * tf;
      }
    diff = std::sqrt( diff );
    if( diff > 2*pnt->GetRadiusInObjectSpace() && diff > 4 )
      {
      std::cout << "*** FAILURE: Local ridge test failed.  Local ridge too far."
        << std::endl;
      std::cout << "   Source = " << xContI << std::endl;
      std::cout << "   Result = " << xRidgeContI << std::endl;
      ++failures;
      continue;
      }

    std::cout << "Local ridge discovery a success!" << std::endl;
    std::cout << std::endl;
    std::cout << "***** Beginning tube extraction ***** " << std::endl;
{'fraction_non_alphanumeric': 0.15997918969890096, 'fraction_numerical': 0.007413669766534434, 'mean_word_length': 3.0736423841059604, 'pattern_counts': {'":': 0, '<': 367, '<?xml version=': 0, '>': 118, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6077805', 'n_tokens_mistral': 5095, 'n_tokens_neox': 4884, 'n_words': 1120}