itk::PSMCommandLineClass< VDimension > Class Template Reference

This class provides a command line tool to run the Particle Shape Modeling. It runs the optimization as well as Procrustes Registration if the option is in the project parameter file. For now, the tool only uses the PSMEntropyModelFilter for the optimization. More...

#include <itkPSMCommandLineClass.h>

Inheritance diagram for itk::PSMCommandLineClass< VDimension >:

Collaboration diagram for itk::PSMCommandLineClass< VDimension >:

Public Types
typedef PSMCommandLineClass	Self
typedef DataObject	Superclass
typedef SmartPointer< Self >	Pointer
typedef SmartPointer< const Self >	ConstPointer
typedef itk::Image< float, VDimension >	ImageType
typedef PSMEntropyModelFilter< typename PSMCommandLineClass::ImageType >	EntropyModelFilterType
typedef PSMProcrustesRegistration< VDimension >	ProcrustesRegistrationType
typedef PSMProjectReader	ProjectReaderType
typedef PSMProject	ProjectType

Public Member Functions
	itkNewMacro (Self)
	itkTypeMacro (PSMCommandLineClass, DataObject)
void	ReadProjectFile (const char *fn)
void	ReadDistanceTransforms (std::string input_path_prefix)
void	ReadCuttingPlanes ()
void	ReadModelPointFiles ()
void	SetPSMFilterParameters ()
void	WriteOutputs (std::string output_path)
void	Run (const char *fname, std::string input_path_prefix, std::string output_path)
	PSMCommandLineClass ()
	itkGetObjectMacro (PSMFilter, EntropyModelFilterType)
	itkGetMacro (ReportInterval, unsigned int)
	itkSetMacro (ReportInterval, unsigned int)

Protected Member Functions
void	IterateCallback (itk::Object *, const itk::EventObject &)

Detailed Description

template<unsigned int VDimension>
class itk::PSMCommandLineClass< VDimension >

This class provides a command line tool to run the Particle Shape Modeling. It runs the optimization as well as Procrustes Registration if the option is in the project parameter file. For now, the tool only uses the PSMEntropyModelFilter for the optimization.

Definition at line 47 of file itkPSMCommandLineClass.h.

Member Typedef Documentation

template<unsigned int VDimension>

typedef PSMEntropyModelFilter<typename PSMCommandLineClass::ImageType> itk::PSMCommandLineClass< VDimension >::EntropyModelFilterType

PSM model optimization filter typedef

Definition at line 67 of file itkPSMCommandLineClass.h.

template<unsigned int VDimension>

typedef itk::Image<float, VDimension> itk::PSMCommandLineClass< VDimension >::ImageType

Input distance transforms image typedef

Definition at line 63 of file itkPSMCommandLineClass.h.

template<unsigned int VDimension>

typedef PSMProcrustesRegistration<VDimension> itk::PSMCommandLineClass< VDimension >::ProcrustesRegistrationType

Procrustes Registration typedef

Definition at line 70 of file itkPSMCommandLineClass.h.

template<unsigned int VDimension>

typedef PSMProjectReader itk::PSMCommandLineClass< VDimension >::ProjectReaderType

Project Reader typedef

Definition at line 73 of file itkPSMCommandLineClass.h.

template<unsigned int VDimension>

typedef PSMProject itk::PSMCommandLineClass< VDimension >::ProjectType

Project typedef

Definition at line 76 of file itkPSMCommandLineClass.h.

template<unsigned int VDimension>

typedef PSMCommandLineClass itk::PSMCommandLineClass< VDimension >::Self

Standard class typedefs.

Definition at line 51 of file itkPSMCommandLineClass.h.

Constructor & Destructor Documentation

template<unsigned int VDimension>

itk::PSMCommandLineClass< VDimension >::PSMCommandLineClass

(

)

Constructor and destructor

Definition at line 33 of file itkPSMCommandLineClass.hxx.

  {
    m_ProcrustesCounter = 0;
    m_ReportInterval = 10;
    m_ProjectReader = ProjectReaderType::New();
    m_Project = ProjectType::New();
    m_ProcrustesRegistration = ProcrustesRegistrationType::New();
    m_PSMFilter = EntropyModelFilterType::New();


  }

Member Function Documentation

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::IterateCallback ( itk::Object *  caller, const itk::EventObject &    )

protected

Callback to run Procrustes Registration on the shapes at the interval specified in the project parameter file or by default.

Definition at line 47 of file itkPSMCommandLineClass.hxx.

  {
    EntropyModelFilterType *o =
      static_cast<EntropyModelFilterType *>(caller);

    // Print every 10 iterations
    if (o->GetNumberOfElapsedIterations() % m_ReportInterval != 0) { return; }

    std::cout << "Iteration # " << o->GetNumberOfElapsedIterations() << std::endl;
    std::cout << " Eigenmode variances: ";
    for (unsigned int i = 0; i < o->GetShapePCAVariances().size(); i++)
    {
      std::cout << o->GetShapePCAVariances()[i] << " ";
    }
    std::cout << std::endl;
    std::cout << " Regularization = " << o->GetRegularizationConstant() << std::endl;

    // Check if optimization is run using scales. Get Procrustes
    // interval for the current scale.
    int interval = 0;
    std::cout << "Optimization Scale " << (o->GetCurrentScale() + 1) << "/"
      << o->GetNumberOfScales() << std::endl;

    if (o->GetNumberOfScales() > 1)
    {
      interval = m_ProcrustesInterval[o->GetCurrentScale()];
      std::cout << "Interval = " << interval << std::endl;
    } else // Use the default interval
    {
      interval = m_ProcrustesRegistration->GetProcrustesInterval();
    }

    // Check if the Procrustes interval is set to a value other than 0
    if (interval > 0)
    {
      this->m_ProcrustesCounter++;

      // If the counter is greater than the interval value, run
      // Procrustes registration
      if (this->m_ProcrustesCounter >= interval)
      {
        // Reset the counter
        this->m_ProcrustesCounter = 0;
        this->m_ProcrustesRegistration->RunRegistration();
        std::cout << "Ran Procrustes Registration" << std::endl;
      }
    }
  }

template<unsigned int VDimension>

itk::PSMCommandLineClass< VDimension >::itkGetMacro	(	ReportInterval	,
		unsigned	int
	)

Get/Set the number of optimization iterations between reporting of the optimized values.

template<unsigned int VDimension>

itk::PSMCommandLineClass< VDimension >::itkGetObjectMacro	(	PSMFilter	,
		EntropyModelFilterType
	)

Returns the particle system used in the shape model computation.

template<unsigned int VDimension>

itk::PSMCommandLineClass< VDimension >::itkNewMacro

(

Self

)

Method for creation through the object factory.

template<unsigned int VDimension>

itk::PSMCommandLineClass< VDimension >::itkTypeMacro	(	PSMCommandLineClass< VDimension >	,
		DataObject
	)

Run-time type information (and related methods).

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::ReadCuttingPlanes

(

)

Reads any cutting plane data from the domain fields in the project file and loads this information into the PSM Filter.

Definition at line 138 of file itkPSMCommandLineClass.hxx.

  {
    if (!m_Project->HasDomains())
    {
      itkExceptionMacro("Domains must be defined in the project file.");
    }

    // Compile the list of domains
    std::vector<std::string> domain_names = m_Project->GetDomainNames();

    for (unsigned int i = 0; i < domain_names.size(); i++)
    {
      if (m_Project->HasDomainCuttingPlanes(domain_names[i]))
      {
        // Cutting planes are returned as a list of vnl_vectors.
        std::vector<vnl_vector_fixed<double, 3> > planes
          = m_Project->GetDomainCuttingPlanes(domain_names[i]);

        // The list should have size that is multiple of 3, since
        // every plane is defined by exactly 3 points.
        if (planes.size() % 3 != 0)
        {
          itkExceptionMacro("Something is wrong with the cutting plane data for domain " + domain_names[i] + ".  Correct number of points was not found.");
        }

        // Traverse the list and pull out sets of 3 points.  Each
        // set of 3 defines a cutting plane.
        vnl_vector_fixed<double, 3> x;
        vnl_vector_fixed<double, 3> y;
        vnl_vector_fixed<double, 3> z;
        std::vector<vnl_vector_fixed<double, 3> >::iterator it =
          planes.begin();

        while (it != planes.end())
        {
          x = *it; it++;
          y = *it; it++;
          z = *it; it++;

          std::cout << "Found cutting plane for domain " << domain_names[i] <<
            ": " << std::endl;
          std::cout << x[0] << " " << x[1] << " " << x[2] << std::endl;
          std::cout << y[0] << " " << y[1] << " " << y[2] << std::endl;
          std::cout << z[0] << " " << z[1] << " " << z[2] << std::endl;

          // Add cutting plane to the PSM filter by domain name

          m_PSMFilter->AddCuttingPlane(x, y, z, domain_names[i]);

          //         int d = m_PSMFilter->GetDomainIndexByName(domain_names[i]);
          //         m_PSMFilter->GetDomain(d)->SetCuttingPlane(x,y,z);

        }


      } // end has cutting planes
    } // end traverse domains

  }

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::ReadDistanceTransforms

(

std::string

input_path_prefix

)

Read the distance transforms that are provided as inputs to the optimzation filter.

Definition at line 98 of file itkPSMCommandLineClass.hxx.

  {
    if (!m_Project->HasDomains())
    {
      throw itk::ExceptionObject("Domains must be defined in the project file.");
    }

    // Compile the list of distance transforms
    std::vector<std::string> domain_names = m_Project->GetDomainNames();
    std::vector<std::string> dt_files;
    for (unsigned int i = 0; i < domain_names.size(); i++)
    {
      std::string s = (m_Project->GetDomainDistanceTransform(domain_names[i]))[0];
      dt_files.push_back(s);


      this->m_PSMFilter->SetDomainName(domain_names[i], i);
    }

    // Read in the distance transforms
    //  const std::vector<std::string> &dt_files = this->m_Project->GetDistanceTransforms();
    std::cout << "Reading distance transforms ..." << std::endl;
    for (unsigned int i = 0; i < dt_files.size(); i++)
    {
      typename itk::ImageFileReader<PSMCommandLineClass::ImageType>::Pointer reader =
        itk::ImageFileReader<PSMCommandLineClass::ImageType>::New();
      reader->SetFileName(input_path_prefix + dt_files[i]);
      //reader->SetFileName(dt_files[i]);
      reader->Update();

      std::cout << domain_names[i] << ": " << dt_files[i] << std::endl;

      this->m_PSMFilter->SetInput(i, reader->GetOutput());

    }
    std::cout << "Done!" << std::endl;
  }

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::ReadModelPointFiles

(

)

Read the input point files, if any.

Definition at line 200 of file itkPSMCommandLineClass.hxx.

  {
    // double value1, value2, value3;
    // unsigned int numOfPoints;
    // typename PSMCommandLineClass::EntropyModelFilterType::PointType pt;

    // // Check if model initialization is supplied. If not, then optimization scales will need to be
    // // used later.
    // if(this->m_Project->HasModel(std::string("initialization")))
    // {
    //   const std::vector<std::string> &pt_files = this->m_Project->GetModel(std::string("initialization"));
    //   std::cout << "Reading the initial model correspondences ..." << std::endl;
    //   for (unsigned int i = 0; i < pt_files.size(); i++)
    //   {
    //     // Read the points for this file and add as a list
    //     typename std::vector<typename PSMCommandLineClass::EntropyModelFilterType::PointType> init_points;
    //     numOfPoints = 0;
    //     // Open the ascii file.
    //     std::ifstream in( (input_path_prefix + pt_files[i]).c_str() );
    //     if ( !in )
    //     {
    //       std::cerr << "Could not open point file for input." << std::endl;
    //       break;
    //     }

    //     // Read all of the points, one point per line.
    //     while (in)
    //     {
    //       in>>value1>>value2>>value3;
    //       pt[0] = value1;
    //       pt[1] = value2;
    //       if(VDimension == 3)
    //       {
    //         pt[2] = value3;
    //       }
    //       init_points.push_back(pt);
    //       numOfPoints++;
    //     }
    //     // this algorithm pushes the last point twice
    //     init_points.pop_back();
    //     std::cout << "Read " << numOfPoints-1 << " points. " << std::endl;
    //     in.close();
    //     this->m_PSMFilter->SetInputCorrespondencePoints(i,init_points);
    //     std::cout << "  " << pt_files[i] << std::endl;
    //   }
    //   std::cout << "Done!" << std::endl;


    // }
    // else
    // {
    //   // Check if scales have been supplied
    //   if(this->m_Project->HasOptimizationAttribute("number_of_scales")
    //      &&
    //      (this->m_Project->GetNumberOfOptimizationScales()) > 0)
    //   {
    //     // Read the input optimization scales from the project file.
    //     this->ReadOptimizationScales();
    //   }
    //   else
    //   {
    //     // REQUIRE THAT WE HAVE SCALES DEFINED.  If there is only a single
    //     // optimization scale, it should be defined as scale 0, but for now
    //     // we are enforcing that scales be defined in the parameter file.

    //     throw itk::ExceptionObject("The number of optimization scales and their parameters were not defined in the parameter file.");

    //     //this->SetDefaultScales();
    //   }
    // }
  }

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::ReadProjectFile

(

const char *

fn

)

Parse the project file and store in the m_Project member variable.

Definition at line 274 of file itkPSMCommandLineClass.hxx.

  {
    // Read the project parameter file
    m_ProjectFileName = fn;
    this->m_ProjectReader->SetFileName(fn);
    this->m_ProjectReader->Update();

    // Store the project parameters  
    this->m_Project = m_ProjectReader->GetOutput();
  }

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::Run	(	const char *	fname,
		std::string	input_path_prefix,
		std::string	output_path
	)

Run the steps of the optimization process

Definition at line 404 of file itkPSMCommandLineClass.hxx.

  {
    // Read the project parameter files.  Also sets the project
    // filename (m_ProjectFileName)
    this->ReadProjectFile(fname);

    // Set the optimization parameters.  This also checks to be sure
    // everything is there.  If not, this method throws an exception.
    this->SetPSMFilterParameters();

    // Read the distance transform files
    this->ReadDistanceTransforms(input_path_prefix);

    // Read any geometric contraints (cutting planes)
    this->ReadCuttingPlanes();

    // Create the callback function to run Procrustes and report
    // progress.
    m_IterateCmd = itk::MemberCommand<PSMCommandLineClass>::New();
    m_IterateCmd->SetCallbackFunction(this, &PSMCommandLineClass::IterateCallback);
    this->m_PSMFilter->AddObserver(itk::IterationEvent(), m_IterateCmd);

    // Set ParticleSystem for ProcrustesRegistration
    this->m_ProcrustesRegistration->SetPSMParticleSystem(m_PSMFilter->GetParticleSystem());

    // Run the optimization filter
    this->m_PSMFilter->Update();
    // TODO: Is this error message really required?
    if (this->m_PSMFilter->GetNumberOfElapsedIterations() >= this->m_PSMFilter->GetMaximumNumberOfIterations())
    {
      std::cerr << "Optimization did not converge based on tolerance criteria." << std::endl;
    }

    std::cout << "Optimization done " << std::endl;

    // Write out the optimized point sets
    this->WriteOutputs(output_path);
  }

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::SetPSMFilterParameters

(

)

Sets the optimization parameters. Assumes that the parameter file has been parsed.

Definition at line 287 of file itkPSMCommandLineClass.hxx.

  {
    // Read the optimization parameters and set their corresponding
    // parameters in the PSM optimization filter.  All of the following
    // parameters are required to be present in the project file.
    std::cout << "Reading optimization parameters ... " << std::endl;
    unsigned int number_of_scales = this->m_Project->GetNumberOfOptimizationScales();
    std::cout << " Optimization scales = " << number_of_scales << std::endl;

    if (number_of_scales > 0) {
      // We need vectors of parameters, one for each scale.
      std::vector<double> regularization_initial(number_of_scales);
      std::vector<double> regularization_final(number_of_scales);
      std::vector<double> regularization_decayspan(number_of_scales);
      std::vector<double> tolerance(number_of_scales);
      std::vector<unsigned int> maximum_iterations(number_of_scales);
      m_ProcrustesInterval.resize(number_of_scales);

      // Read parameters for each scale
      for (unsigned int i = 0; i < number_of_scales; i++)  {
        std::cout << "Optimization parameters for scale " << i << ": " << std::endl;
        if (this->m_Project->HasOptimizationAttribute("regularization_initial", i))
        {
          regularization_initial[i] = this->m_Project->GetOptimizationAttribute("regularization_initial", i);
        } else { throw itk::ExceptionObject("Missing parameter: regularization_initial"); }
        std::cout << "    regularization_initial = " << regularization_initial[i] << std::endl;

        if (this->m_Project->HasOptimizationAttribute("regularization_final", i))
        {
          regularization_final[i] = this->m_Project->GetOptimizationAttribute("regularization_final", i);
        } else { throw itk::ExceptionObject("Missing parameter: regularization_final"); }
        std::cout << "      regularization_final = " << regularization_final[i] << std::endl;

        if (this->m_Project->HasOptimizationAttribute("regularization_decayspan", i))
        {
          regularization_decayspan[i]
            = this->m_Project->GetOptimizationAttribute("regularization_decayspan", i);
        } else { throw itk::ExceptionObject("Missing parameter: regularization_decayspan"); }
        std::cout << "  regularization_decayspan = " << regularization_decayspan[i] << std::endl;

        if (this->m_Project->HasOptimizationAttribute("tolerance", i))
        {
          tolerance[i] = this->m_Project->GetOptimizationAttribute("tolerance", i);
        } else { throw itk::ExceptionObject("Missing parameter: tolerance"); }
        std::cout << "                 tolerance = " << tolerance[i] << std::endl;

        if (this->m_Project->HasOptimizationAttribute("maximum_iterations", i))
        {
          maximum_iterations[i]
            = static_cast<unsigned int>(this->m_Project->GetOptimizationAttribute("maximum_iterations", i));
        } else { throw itk::ExceptionObject("Missing parameter: maximum_iterations"); }
        std::cout << "        maximum_iterations = " << maximum_iterations[i] << std::endl;

        //       if(this->m_Project->HasProcrustes() == true)
        //   {
        if (this->m_Project->HasOptimizationAttribute("procrustes_interval", i))
        {
          m_ProcrustesInterval[i]
            = static_cast<unsigned int>(this->m_Project->GetOptimizationAttribute("procrustes_interval", i));
        } else { throw itk::ExceptionObject("Missing parameter: procrustes_interval"); }
        std::cout << "        procrustes_interval = " << m_ProcrustesInterval[i] << std::endl;
        //}      
      }

      // Now set the parameters in the filter
      this->m_PSMFilter->SetNumberOfScales(number_of_scales);
      this->m_PSMFilter->SetRegularizationInitial(regularization_initial);
      this->m_PSMFilter->SetRegularizationFinal(regularization_final);
      this->m_PSMFilter->SetRegularizationDecaySpan(regularization_decayspan);
      this->m_PSMFilter->SetTolerance(tolerance);
      this->m_PSMFilter->SetMaximumNumberOfIterations(maximum_iterations);
    } else
    {
      throw itk::ExceptionObject("Project must define at least one optimization scale");
    }
  }

template<unsigned int VDimension>

void itk::PSMCommandLineClass< VDimension >::WriteOutputs

(

std::string

output_path

)

Write out the optimized point sets to user specified files

Definition at line 366 of file itkPSMCommandLineClass.hxx.

  {
    const std::vector<std::string> &out_files = this->m_Project->GetModel(std::string("optimized"));
    if (out_files.size() != this->m_PSMFilter->GetParticleSystem()->GetNumberOfDomains())
    {
      std::cerr << "Number of output files does not match the number of particle domains (inputs)." << std::endl;
    } else
    {
      for (unsigned int d = 0; d < this->m_PSMFilter->GetParticleSystem()->GetNumberOfDomains(); d++)
      {
        // Open the output file.
        std::string fname = output_path + out_files[d];
        std::ofstream out(fname.c_str());
        if (!out)
        {
          std::cerr << "Could not open point file for output: " << std::endl;
          break;
        } else
        {
          for (unsigned int j = 0; j < this->m_PSMFilter->GetParticleSystem()->GetNumberOfParticles(d); j++)
          {
            for (unsigned int i = 0; i < VDimension; i++)
            {
              out << this->m_PSMFilter->GetParticleSystem()->GetPosition(j, d)[i] << " ";
            }
            if (VDimension == 2)
            {
              out << "0.0"; // printed for display in ShapeWorksViewer2
            }
            out << std::endl;
          }
        }
      } // end for number of domains
    }
  }

---

The documentation for this class was generated from the following files:

• C:/Users/Brig/Documents/ShapeWorksStudio/src/ParticleShapeworks/Examples/itkPSMCommandLineClass.h
• C:/Users/Brig/Documents/ShapeWorksStudio/src/ParticleShapeworks/Examples/itkPSMCommandLineClass.hxx