C++ GetPot::follow方法代码示例

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }
  const double a = cl.follow (0.0, "-a");
  const double b = cl.follow (1.0, "-b");
  const unsigned int nnodes = cl.follow (100, 2, "-n", "--nnodes");
  const std::string diffusion = cl.follow ("1.0", 2, "-d", "--diffusion");
  const std::string forcing = cl.follow ("1.0", 2, "-f", "--forcing");
  
  std::unique_ptr<mesh<double>> m (new mesh<double> (a, b, nnodes));
  fem_1d<double> problem (std::move (m));

  coeff<double> a_coeff (diffusion);
  problem.set_diffusion_coefficient (a_coeff);

  coeff<double> f_coeff (forcing);
  problem.set_source_coefficient (f_coeff);

  problem.assemble ();
  
  problem.set_dirichlet (fem_1d<double>::left_boundary, 0.0);
  problem.set_dirichlet (fem_1d<double>::right_boundary, 0.0);

  problem.solve ();
  
  for (unsigned int ii = 0; ii < nnodes; ++ii)
    std::cout << problem.m->nodes[ii] << " "
              << problem.result ()(ii, 0)
              << std:: endl;
      
  return 0;
};

int main (int argc, char *argv[])
	{
	const unsigned int ImageDimension = 2;
	    GetPot cl (argc, const_cast<char**>(argv));
	    if( cl.size() == 1 || cl.search (2,"--help","-h") )
	    {
	        std::cout << "Not Enough Arguments" << std::endl;
	        std::cout << "Generate the Gradient Table" << std::endl;
	        std::cout << "Usage:  return -1" << std::endl;
	    }

	    const string image_n = cl.follow("NoFile",1, "-i");
	    const string out_n   = cl.follow("NoFile",1, "-o");

	    typedef itk::DiffusionTensor3D<float> DiffusionTensorType;
	    typedef itk::Image<DiffusionTensorType, 3> ImageType;
	    typedef itk::ImageFileReader<ImageType> ReaderType;
	    ReaderType::Pointer reader = ReaderType::New();

	    reader->SetFileName(image_n);
	    reader->Update();

	    ImageType::Pointer image = reader->GetOutput();

	    typedef itk::ImageRegionIterator<ImageType> TensorIterator;
	    TensorIterator itImg(image, image->GetLargestPossibleRegion());

	    std::ofstream  file;
	    file.open(out_n);

	    for(itImg.GoToBegin(); !itImg.IsAtEnd(); ++itImg)
	    {
	    	file << itImg.Get() << std::endl;
	    }

	    file.close();



		return 0;
	}

int main (int argc, char **argv)
{

    GetPot cl (argc, argv);

    if (cl.search (2, "-h", "--help"))
    {
        std::cerr << help_text << std::endl;
        return 0;
    }
    const double a = cl.follow (0.0, "-a");
    const double b = cl.follow (1.0, "-b");
    const unsigned int nnodes = cl.follow (100, 2, "-n", "-nnodes");

    mesh m (a, b, nnodes);
    Eigen::MatrixXd A(nnodes, nnodes);

    Eigen::Matrix2d mloc;
    mloc << 0, 0, 0, 0;
    for (unsigned int iel = 0; iel < m.nels; ++iel)
    {

        mloc << 0, 0, 0, 0;
        for (unsigned int inode = 0; inode < 2; ++inode)
        {
            double igrad = (inode == 0 ? 1.0 / m.h : -1.0 / m.h);
            for (unsigned int jnode = 0; jnode < 2; ++jnode)
            {
                double jgrad =  (jnode == 0 ? 1.0 / m.h : -1.0 / m.h);
                mloc(inode,jnode) = igrad * jgrad * m.h;
                A(m.elements[iel][inode],m.elements[iel][jnode]) +=
                    mloc(inode,jnode);
            }
        }
    }

    Eigen::VectorXd f(nnodes);
    Eigen::Vector2d vloc;

    for (unsigned int iel = 0; iel < m.nels; ++iel)
    {
        vloc << 0, 0;

        for (unsigned int inode = 0; inode < 2; ++inode)
        {
            vloc(inode) = m.h / 2.0;
            f(m.elements[iel][inode]) += vloc(inode);
        }
    }

    f(0) = 0;
    f(nnodes - 1) = 0;

    A(0,0) = 1.0;
    A(nnodes-1,nnodes-1) = 1.0;
    for (unsigned int ii = 1; ii < nnodes; ++ii)
    {
        A(0, ii) = 0.0;
        A(nnodes-1, nnodes-1-ii) = 0.0;
    }


    Eigen::VectorXd uh = A.partialPivLu ().solve (f);

    for (unsigned int ii = 0; ii < nnodes; ++ii)
        std::cout << uh(ii, 0) << std:: endl;

    return 0;
};

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }
  const double a = cl.follow (0.0, "-a");
  const double b = cl.follow (1.0, "-b");
  const unsigned int nnodes = cl.follow (100, 2, "-n", "-nnodes");

  mesh m (a, b, nnodes);
  matrix A(nnodes);

  matrix mloc(2);
  for (unsigned int iel = 0; iel < m.nels; ++iel)
    {

      std::fill (mloc.get_data (),
                 mloc.get_data () + 4,
                 0.0);
      
      for (unsigned int inode = 0; inode < 2; ++inode)
        {
          double igrad = (inode == 0 ? 1.0 / m.h : -1.0 / m.h);
          for (unsigned int jnode = 0; jnode < 2; ++jnode)
            {
              double jgrad =  (jnode == 0 ? 1.0 / m.h : -1.0 / m.h);
              mloc(inode,jnode) = igrad * jgrad * m.h;
              A(m.elements[iel][inode],m.elements[iel][jnode]) += 
                mloc(inode,jnode);
            }
        }
    }

  matrix f(nnodes, 1);
  matrix vloc(2, 1);
  
  for (unsigned int iel = 0; iel < m.nels; ++iel)
    {
      std::fill (mloc.get_data (),
                 mloc.get_data () + 2,
                 0.0);
      
      for (unsigned int inode = 0; inode < 2; ++inode)
        {
          vloc(inode, 0) = m.h / 2.0;
          f(m.elements[iel][inode], 0) += vloc(inode, 0);
        }
    }

  f(0, 0) = 0;
  f(nnodes - 1, 0) = 0;

  A(0,0) = 1.0;
  A(nnodes-1,nnodes-1) = 1.0;
  for (unsigned int ii = 1; ii < nnodes; ++ii)
    {
      A(0, ii) = 0.0;
      A(nnodes-1, nnodes-1-ii) = 0.0;
    }

  
  matrix uh(f);
  A.solve (uh);

  for (unsigned int ii = 0; ii < nnodes; ++ii)
    std::cout << uh(ii, 0) << std:: endl;
      
  return 0;
};

int main (int argc, char** argv)
{
#ifdef LIFEV_HAS_HDF5
#ifdef HAVE_MPI

    typedef RegionMesh<LinearTetra> mesh_Type;

    MPI_Init (&argc, &argv);
    boost::shared_ptr<Epetra_Comm> comm (new Epetra_MpiComm (MPI_COMM_WORLD) );

    if (comm->NumProc() != 1)
    {
        std::cout << "This test needs to be run "
                  << "with a single process. Aborting."
                  << std::endl;
        return (EXIT_FAILURE);
    }

    GetPot commandLine (argc, argv);
    string dataFileName = commandLine.follow ("data", 2, "-f", "--file");
    GetPot dataFile (dataFileName);

    const UInt numElements (dataFile ("mesh/nelements", 10) );
    const UInt numParts (dataFile ("test/num_parts", 4) );
    const std::string partsFileName (dataFile ("test/hdf5_file_name", "cube.h5") );
    const std::string ioClass (dataFile ("test/io_class", "new") );

    std::cout << "Number of elements in mesh: " << numElements << std::endl;
    std::cout << "Number of parts: " << numParts << std::endl;
    std::cout << "Name of HDF5 container: " << partsFileName << std::endl;

    boost::shared_ptr<mesh_Type> fullMeshPtr (new mesh_Type ( comm ) );
    regularMesh3D (*fullMeshPtr, 1, numElements, numElements, numElements,
                   false, 2.0, 2.0, 2.0, -1.0, -1.0, -1.0);

    MeshPartitioner<mesh_Type> meshPart;
    meshPart.setup (numParts, comm);

    meshPart.attachUnpartitionedMesh (fullMeshPtr);
    meshPart.doPartitionGraph();
    meshPart.doPartitionMesh();

    // Release the original mesh from the MeshPartitioner object and
    // delete the RegionMesh object
    meshPart.releaseUnpartitionedMesh();
    fullMeshPtr.reset();

    // Write mesh parts to HDF5 container
    if (! ioClass.compare ("old") )
    {
        ExporterHDF5Mesh3D<mesh_Type> HDF5Output (dataFile,
                                                  meshPart.meshPartition(),
                                                  partsFileName,
                                                  comm->MyPID() );
        HDF5Output.addPartitionGraph (meshPart.elementDomains(), comm);
        HDF5Output.addMeshPartitionAll (meshPart.meshPartitions(), comm);
        HDF5Output.postProcess (0);
        HDF5Output.closeFile();
    }
    else
    {
        boost::shared_ptr<Epetra_MpiComm> mpiComm =
            boost::dynamic_pointer_cast<Epetra_MpiComm> (comm);
        PartitionIO<mesh_Type> partitionIO (partsFileName, mpiComm);
        partitionIO.write (meshPart.meshPartitions() );
    }

    MPI_Finalize();

#else
    std::cout << "This test needs MPI to run. Aborting." << std::endl;
    return (EXIT_FAILURE);
#endif /* HAVE_MPI */
#else
    std::cout << "This test needs HDF5 to run. Aborting." << std::endl;
    return (EXIT_FAILURE);
#endif /* LIFEV_HAS_HDF5 */

    return (EXIT_SUCCESS);
}

void Server::parseConfig(GetPot &config)
{
  // this variable is needed for setting up the partial loading part when the features
  // are loaded for the very first time at fire start up
  string partialLoadingString="empty";

  if(config.search(2,"-F","--filter"))
  {
    string unparsed = config.follow("empty",2,"-F","--filter");
    if (parseFilter(unparsed.c_str()))
    {
      retriever_.setFilterApply(true);
      DBG(10) << "Using filter " << unparsed << " for retrieval" << endl;
    }
    else
    {
      retriever_.setFilterApply(false);
      retriever_.clearFilter();
      DBG(10) << "Error in filter String " << unparsed << endl;
      DBG(10) << "Aborting filtering and using normal retrieval settings" << endl;
    }
  }

  if(config.search(2,"-q","--queryCombiner")) {
    retriever_.setQueryCombiner(config.follow("adding:POS=1.0:NEG=0.833",2,"-q","--queryCombiner"));
  }

  if(config.search("--reRanker")) {
    retriever_.setReranking(config.follow("cluster:CONS=100:RR=20:CLUSTERS=5","--reRanker"));
  } 
  if(config.search(2,"-U","--defdontload"))
  {
    partialLoadingString="default";
    retriever_.setPartialLoadingApply(true);
  }

  if(config.search(2,"-u","--dontload"))
  {
    string unparsed=config.follow("empty",2,"-u","--dontload");
    if(unparsed!="empty")
    {
      // the parsing is not possible here, because the amount of different features
      // is not known yet
      // but a sanity check is possible
      const char* str = unparsed.c_str();
      bool errformat = false;
      // note it is not checked if there are two or more : present in a row
      while(*str != '\0' && !errformat)
      {
        if(!isdigit(*str) && *str!=':')
        {
          errformat=true;
          DBG(10) << "Error in format string " << unparsed << " following -u/--dontload option" << endl;
          DBG(10) << "format string must start with a digit; aborting partial feature loading" << endl;
        }
        str++;
      }
      if(!errformat && retriever_.setPartialLoadingApply(true))
      {
        partialLoadingString=unparsed;
      }
      else
      {
        ERR << "-u/--dontload was used without -F/--filter option or there are illegal characters in the format string" << endl;
        ERR << "ignoring -u/--dontload option, loading all features" << endl;
      }
    }
    else
    {
      retriever_.setPartialLoadingApply(false);
      ERR << "Error after -u/--dontload option: No featureindexsequence given" << endl;
      //DBG(10) << "Aborting partial loading for filtered retrieval and loading all features" << endl;
    }
  }

  if(config.search(2,"-f","--filelist"))
  {
    string filelistname=config.follow("list.txt",2,"-f","--filelist");
    string result=retriever_.filelist(filelistname,partialLoadingString);
    if(result=="filelist FAILURE")
    {
      ERR << "Error in filelist. Please correct it and start again." << endl;
      exit(1);
    }
    DBG(10) << result << endl;
  }
  if(config.search(2,"-R","--relevancefile"))
  {
    relevanceFile_=config.follow("",2,"-R","--relevancefile");
  }

  retriever_.setScoring("linear");

  
  log_=LogFile(config.follow("",2,"-l","--logfile"));

  if(config.search(2,"-D","--defaultdists"))
  {
    for(uint i=0;i<retriever_.numberOfSuffices();++i)
    {
//.........这里部分代码省略.........

int
main ( int argc, char** argv )
{
#ifdef LIFEV_HAS_HDF5
#ifdef HAVE_MPI

    MPI_Init (&argc, &argv);
    boost::shared_ptr<Epetra_MpiComm> comm (new Epetra_MpiComm (MPI_COMM_WORLD) );

    const bool verbose (comm->MyPID() == 0);

    // Read first the data needed

    if (verbose)
    {
        std::cout << " -- Reading the data ... " << std::flush;
    }
    //    GetPot dataFile ( "data" );
    if (verbose)
    {
        std::cout << " done ! " << std::endl;
    }

    GetPot cl (argc, argv);
    // partitionerType should be MeshPartitioner, MeshPartitionTool_ParMETIS or
    // MeshPartitionTool_Zoltan
    const std::string partitionerType = cl.follow ("MeshPartitioner",
                                                   "--partitioner-type");
    std::string partsFile;
    partsFile.reserve (50);
    partsFile += "cube_";
    partsFile += partitionerType;
    partsFile += ".h5";

    boost::shared_ptr<mesh_Type> mesh;
    {
        PartitionIO<RegionMesh<LinearTetra> > partitionIO (partsFile, comm);
        partitionIO.read (mesh);
    }

    // Build the FESpaces

    if (verbose)
    {
        std::cout << " -- Building FESpaces ... " << std::flush;
    }
    std::string uOrder ("P1");
    std::string bOrder ("P1");
    boost::shared_ptr<FESpace<mesh_Type, MapEpetra> >
    uFESpace (new FESpace<mesh_Type, MapEpetra> (mesh, uOrder, 1, comm) );
    boost::shared_ptr<FESpace<mesh_Type, MapEpetra> >
    betaFESpace (new FESpace<mesh_Type, MapEpetra> (mesh, bOrder, 3, comm) );
    if (verbose)
    {
        std::cout << " done ! " << std::endl;
    }
    if (verbose) std::cout << " ---> Dofs: "
                               << uFESpace->dof().numTotalDof() << std::endl;

    // Build the assembler and the matrices

    if (verbose)
    {
        std::cout << " -- Building assembler ... " << std::flush;
    }
    ADRAssembler<mesh_Type, matrix_Type, vector_Type> adrAssembler;
    if (verbose)
    {
        std::cout << " done! " << std::endl;
    }

    if (verbose)
    {
        std::cout << " -- Setting up assembler ... " << std::flush;
    }
    adrAssembler.setup (uFESpace, betaFESpace);
    if (verbose)
    {
        std::cout << " done! " << std::endl;
    }

    if (verbose)
    {
        std::cout << " -- Defining the matrix ... " << std::flush;
    }
    boost::shared_ptr<matrix_Type>
    systemMatrix (new matrix_Type (uFESpace->map() ) );
    *systemMatrix *= 0.0;
    if (verbose)
    {
        std::cout << " done! " << std::endl;
    }

    // Perform the assembly of the matrix

    if (verbose)
    {
        std::cout << " -- Adding the diffusion ... " << std::flush;
    }
    adrAssembler.addDiffusion (systemMatrix, 1);
//.........这里部分代码省略.........

int main (int argc, char *argv[])
{

    GetPot cl (argc, const_cast<char**>(argv));
    if( cl.size() == 1 || cl.search (2,"--help","-h") )
    {
        std::cout << "Not Enough Arguments" << std::endl;
        std::cout << "Scales the tensors with a scalar factor" << std::endl;
        std::cout << "Usage: -trueB0 <true B0> -m <MaskImage> -true <True Tensors> -f <flag for extended gradient> -t <initial tensor estimate> -g <gradient> -o <Output File> -s <Sigma> -nm <Noise Model> -Sim <intelligent COnvergence>" << std::endl;
      return -1;
    }

    const string image_n =cl.follow("NoFile", 1, "-i");
    const string out_n = cl.follow("NoFile", 1, "-o");

    const int idx = cl.follow(0, 1, "-ix");
    const int idy = cl.follow(0, 1, "-iy");
    const int idz = cl.follow(0, 1, "-iz");

    const int sx = cl.follow(0, 1, "-sx");
    const int sy = cl.follow(0, 1, "-sy");
    const int sz = cl.follow(0, 1, "-sz");



    typedef itk::Image< float, 3> ImageType;
    typedef itk::ImageFileReader<ImageType> ReaderType;
    ReaderType::Pointer reader = ReaderType::New();

    reader->SetFileName(image_n);
    reader->Update();

    ImageType::Pointer image = reader->GetOutput();

    ImageType::SizeType size;
    ImageType::IndexType id;

    id[0] = idx; id[1] = idy; id[2] = idz;

    size[0] = sx; size[1] = sy; size[2] =sz;

    ImageType::RegionType region(id, size);

    typedef itk::ExtractImageFilter<ImageType, ImageType> FilterType;
    FilterType::Pointer filter = FilterType::New();

    filter->SetExtractionRegion(region);
    filter->SetInput(image);
    filter->SetDirectionCollapseToIdentity();

    filter->Update();

    ImageType::Pointer exImage = filter ->GetOutput();

    typedef itk::ImageFileWriter<ImageType> WriterType;
    WriterType::Pointer writer = WriterType::New();

    writer->SetFileName(out_n);
    writer->SetInput(exImage);
    writer->Update();






    return 0;
}

int main(int argc, char *argv[])
{
    GetPot cl (argc, const_cast<char**>(argv));
    if( cl.size() == 1 || cl.search (2,"--help","-h") )
	    {
	        std::cout << "Not Enough Arguments" << std::endl;
	        std::cout << "Generate the Gradient Table" << std::endl;
	        std::cout << "Usage:  return -1" << std::endl;
	    }

    // One idea is to apply the rotation matrix
    const string Linear_trans_n = cl.follow("NoFile",1,"-t");
    const string ref_Image_n = cl.follow("NoFile",1, "-r");
    const string out_n = cl.follow("NoFile",1,"-o");
    const string B0_n =cl.follow("NoFile",1,"-B0");
    const string T1_n =cl.follow("NoFile",1,"-T1");

    typedef itk::Image<float, 3> ImageType;
    typedef itk::ImageFileReader<ImageType> ImageReaderType;
    ImageReaderType::Pointer imageReaderB0 = ImageReaderType::New();
    imageReaderB0->SetFileName(B0_n);
    imageReaderB0->Update();

    std::cout << imageReaderB0->GetOutput()->GetDirection() << std::endl;


    ImageReaderType::Pointer imageReaderT1 = ImageReaderType::New();

    imageReaderT1->SetFileName(B0_n);
        imageReaderT1->Update();


        std::cout << imageReaderT1->GetOutput()->GetDirection() << std::endl;


//    typedef itk::TransformFileReader TransformFileReaderType;
//    typedef TransformFileReaderType::TransformListType TransformListType;
//    typedef itk::TransformBase TransformBaseType;
//    typedef itk::AffineTransform<double, 3> AffineTransformType;
//
//    typedef itk::Image<float, 3> ImageType;
//    typedef itk::ImageFileReader<ImageType> ImageReaderType;
//    ImageReaderType::Pointer imageReader = ImageReaderType::New();
//
//    imageReader->SetFileName(ref_Image_n);
//    imageReader->Update();
//    ImageType::Pointer refImage = imageReader->GetOutput();
//
//
//
//	TransformFileReaderType::Pointer readerTransform = TransformFileReaderType::New();
//	readerTransform->SetFileName(Linear_trans_n);
//	readerTransform -> Update();
//	TransformListType *list = readerTransform->GetTransformList();
//	TransformBaseType * transform = list->front().GetPointer();
//	TransformBaseType::ParametersType parameters = transform->GetParameters();
//	AffineTransformType::Pointer transform_fwd = AffineTransformType::New();
//	transform_fwd->SetParameters(parameters);
//
//	std::cout << transform_fwd->GetParameters() << std::endl;

//    typedef itk::Vector< float, 3 >          VectorPixelType;
//    typedef itk::Image< VectorPixelType, 3 > DisplacementFieldImageType;
//
//
//    typedef itk::TransformToDisplacementFieldFilter<DisplacementFieldImageType, double> DisplacementFieldGeneratorType;
//    DisplacementFieldGeneratorType::Pointer dispfieldGenerator = DisplacementFieldGeneratorType::New();
//
//    dispfieldGenerator->UseReferenceImageOn();
//    dispfieldGenerator->SetReferenceImage( refImage );
//    dispfieldGenerator->SetTransform( transform_fwd );
//    dispfieldGenerator->Update();
//    DisplacementFieldImageType::Pointer dispField = dispfieldGenerator->GetOutput();
//
//    typedef itk::ImageFileWriter<DisplacementFieldImageType> DispFieldWriterType;
//    DispFieldWriterType::Pointer writer = DispFieldWriterType::New();
//    writer->SetFileName(out_n);
//    writer->SetInput(dispField);
//    writer->Update();

    return 0;
}

int main (int argc, char *argv[])
{

    GetPot cl (argc, const_cast<char**>(argv));
    if( cl.size() == 1 || cl.search (2,"--help","-h") )
    {
        std::cout << "Not Enough Arguments" << std::endl;
        std::cout << "Scales the tensors with a scalar factor" << std::endl;
        std::cout << "Usage: -trueB0 <true B0> -m <MaskImage> -true <True Tensors> -f <flag for extended gradient> -t <initial tensor estimate> -g <gradient> -o <Output File> -s <Sigma> -nm <Noise Model> -Sim <intelligent COnvergence>" << std::endl;
      return -1;
    }




     const string file_g_n = cl.follow("NoFile",1, "-g");
     const string fileIn  = cl.follow("NoFile",1,"-i");
     const string fileIn_HR  = cl.follow("NoFile",1,"-iHR");
     const string B0_n     =  cl.follow("NoFile", 1, "-B0");
     const string mask_LR_n   = cl.follow("NoFile",1, "-mLR");
     const int numOfIter   = cl.follow(1,1, "-n");
     const float kappa_L   = cl.follow(0.05, 1, "-k");
     const float lambda_L  = cl.follow(0.25, 1, "-lamb_L");
     const string B0Image_HR_n = cl.follow("NoFile", 1, "-B0HR");
     const string dispField_n = cl.follow("NoFile",1,"-d");
     const string T1Image_n   = cl.follow("NoFile",1,"-T1");
     const string mask_HR_n = cl.follow("NoFile", 1, "-mHR");
    // Usual Typedefs
    typedef float RealType;
    const int ImageDim  =3;

    typedef itk::Image<RealType, ImageDim> ScalarImageType;
    typedef itk::Vector<double, ImageDim> VectorType;
    typedef itk::Image<VectorType, ImageDim> DeformationFieldType;
    typedef itk::Image<VectorType, ImageDim> VectorImageType;

    typedef itk::ImageFileReader<ScalarImageType> ScalarImageReaderType;
    typedef itk::ImageFileWriter<ScalarImageType> ScalarImageWriterType;

    //Read T1 image
	ScalarImageReaderType::Pointer scalarReader = ScalarImageReaderType::New();
	scalarReader->SetFileName(T1Image_n.c_str());
	scalarReader->Update();
	
	ScalarImageType::Pointer T1_image = scalarReader->GetOutput();

    // Read LR ImageList
	typedef std::vector<ScalarImageType::Pointer> ImageListType;
	ImageListType DWIList;
	
	std::ifstream file(fileIn.c_str());
        int numOfImages = 0;
        file >> numOfImages;	
    
	for (int i=0; i < numOfImages  ; i++) // change of numOfImages
         {
             char filename[256];
             file >> filename;
             ScalarImageReaderType::Pointer myReader=ScalarImageReaderType::New();
             myReader->SetFileName(filename);
             std::cout << "Reading.." << filename << std::endl; // add a try catch block
             myReader->Update();
             DWIList.push_back( myReader->GetOutput() ); //using push back to create a stack of diffusion images
         }

	// Read deformation field
	typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
	DeformationFieldReaderType::Pointer deformationFieldReader = DeformationFieldReaderType::New();

	deformationFieldReader->SetFileName(dispField_n.c_str());
	deformationFieldReader->Update();
	
	DeformationFieldType::Pointer defField = deformationFieldReader->GetOutput();
	

	// Read Mask Image Spatial
	typedef itk::ImageMaskSpatialObject<ImageDim> MaskSpatialObjectType;
	typedef MaskSpatialObjectType::ImageType MaskSpatialImageType;
	typedef itk::ImageFileReader<MaskSpatialImageType> MaskSpatialImageReaderType;	

	MaskSpatialImageReaderType::Pointer spatialReader = MaskSpatialImageReaderType::New();
	spatialReader->SetFileName(mask_LR_n.c_str());
	spatialReader->Update();
	
	MaskSpatialImageType::Pointer maskSpatialImage_LR = spatialReader->GetOutput();

	//Read Mask Image Normal
	ScalarImageReaderType::Pointer maskImageReader = ScalarImageReaderType::New();
	maskImageReader->SetFileName(mask_LR_n.c_str());
	maskImageReader->Update();
	
	ScalarImageType::Pointer maskImage_LR = maskImageReader->GetOutput();

		
	// Resample diffusion Images	
/*	typedef itk::WarpImageFilter<ScalarImageType, ScalarImageType, DeformationFieldType> WarpImageFilterType;
	typedef itk::ImageFileWriter<ScalarImageType> ScalarImageWriterType;
	
	

//.........这里部分代码省略.........

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }

  const double a =
    cl.follow (0.0, "-a");

  const double b =
    cl.follow (1.0, "-b");

  const unsigned int nnodes =
    cl.follow (100, 2, "-n", "--nnodes");

  const std::string diffusion =
    cl.follow ("1.0", 2, "-d", "--diffusion");

  const std::string forcing =
    cl.follow ("1.0", 2, "-f", "--forcing");
  
  coeff f_coeff (forcing);
  coeff a_coeff (diffusion);

  const std::string quadrature =
    cl.follow ("trapezoidal.so",
               2, "-q", "--quadrature-rule");

  std::function <void ()> r_r;
  
  void * handle = dlopen (quadrature.c_str (), RTLD_NOW);
  if (! handle)
    {
      std::cerr << "fem1d: cannot load dynamic object!"
                << std::endl;
      std::cerr << dlerror ()
                << std::endl;
      return (-1);
    }

  void * sym = dlsym (handle, "register_rules");
  if (! sym)
    {
      std::cerr << "fem1d: cannot load symbol!"
                << std::endl;
      std::cerr << dlerror ()
                << std::endl;
      return (-1);
    }


  r_r = reinterpret_cast <void (*) ()> (sym);
  r_r ();

  auto & the_factory = quadrature_factory::instance (); 
  auto integrate = the_factory.create ("trapezoidal");
  if (! integrate)
    {
      std::cerr << "rule name unknown" << std::endl;
      return (-1);
    }
    
  mesh m (a, b, nnodes);
  Eigen::SparseMatrix<double> A(nnodes, nnodes);
  
  Eigen::Matrix2d mloc;
  mloc << 0, 0, 0, 0;
  for (unsigned int iel = 0; iel < m.nels; ++iel)
    {

      mloc << 0, 0, 0, 0;      
      for (unsigned int inode = 0; inode < 2; ++inode)
        {

          auto igrad = [=, &m] (double x) -> double
            {
              return basisfungrad
              (x, m.nodes[m.elements[iel][0]],
               m.nodes[m.elements[iel][1]],
               inode == 0 ? 1.0 : 0.0, 
               inode == 1 ? 1.0 : 0.0);
            };

          
          for (unsigned int jnode = 0; jnode < 2; ++jnode)
            {

              auto jgrad = [=, &m] (double x) -> double
                {
                  return basisfungrad
                  (x, m.nodes[m.elements[iel][0]],
                   m.nodes[m.elements[iel][1]],
                   jnode == 0 ? 1.0 : 0.0, 
                   jnode == 1 ? 1.0 : 0.0);
                };
//.........这里部分代码省略.........

int main (int argc, char *argv[])
	{
	const unsigned int ImageDimension = 2;
	    GetPot cl (argc, const_cast<char**>(argv));
	    if( cl.size() == 1 || cl.search (2,"--help","-h") )
	    {
	        std::cout << "Not Enough Arguments" << std::endl;
	        std::cout << "Generate the Gradient Table" << std::endl;
	        std::cout << "Usage:  return -1" << std::endl;
	    }

	   const string image_n = cl.follow("NoFile",1, "-i");
	   const string mask_n = cl.follow("NoFile", 1, "-m");	 
	   const string out_n   = cl.follow("NoFile",1, "-o");

	    typedef itk::DiffusionTensor3D<float> DiffusionTensorType;
	    typedef itk::Image<DiffusionTensorType, 3> TensorImageType;
	    typedef itk::ImageFileReader<TensorImageType> TensorReaderType;
	    TensorReaderType::Pointer reader = TensorReaderType::New();

	    reader->SetFileName(image_n.c_str());
	    reader->Update();

	    TensorImageType::Pointer image = reader->GetOutput();


	    typedef itk::Image<float, 3> ScalarImageType;
	    typedef itk::ImageFileReader<ScalarImageType> ScalarReaderType;

	    ScalarReaderType::Pointer scalarReader = ScalarReaderType::New();
	    scalarReader->SetFileName(mask_n.c_str());
	    scalarReader->Update();
	    ScalarImageType::Pointer maskImage = scalarReader->GetOutput();
	
	    typedef itk::ImageRegionIterator<TensorImageType> TensorIterator;

	    typedef itk::ImageRegionIterator<ScalarImageType> ScalarIterator;
	    ScalarIterator itMask(maskImage, maskImage->GetLargestPossibleRegion());
	

	    TensorUtilities  utilsTensor;
	    TensorImageType::Pointer logTensorImage = utilsTensor.LogTensorImageFilter(image, maskImage);


	    typedef itk::ImageFileWriter<TensorImageType> TensorImageWriterType;
	    TensorImageWriterType::Pointer tensorImageWriter = TensorImageWriterType::New();

            tensorImageWriter->SetFileName("LogTensorImage_stupid.nii.gz");
	    tensorImageWriter->SetInput(logTensorImage);
	    tensorImageWriter->Update();
	
	    std::ofstream  file;
	    file.open(out_n);

	 	ScalarImageType::Pointer TraceImage = ScalarImageType::New();
	    CopyImage	cpImage;
	    cpImage.CopyScalarImage(maskImage, TraceImage);

	    ScalarIterator itTr(TraceImage, TraceImage->GetLargestPossibleRegion());



	    TensorIterator itImg(logTensorImage, logTensorImage->GetLargestPossibleRegion());
	    for(itImg.GoToBegin(), itMask.GoToBegin(), itTr.GoToBegin(); !itTr.IsAtEnd(), !itImg.IsAtEnd(), !itMask.IsAtEnd(); ++itTr, ++itImg, ++itMask)
	    {
	    		file << itImg.Get().GetTrace() << std::endl;
	    		itTr.Set(itImg.Get().GetTrace()) ;
		

	   }


		typedef itk::ImageFileWriter<ScalarImageType> WriterType;
		WriterType::Pointer writer = WriterType::New();
	
		writer->SetFileName("LogTense_Trace.nii.gz");
		writer->SetInput(TraceImage);
		writer->Update();


	    file.close();


		return 0;
	}

int main (int argc, char **argv)
{

  GetPot cl (argc, argv);

  if (cl.search (2, "-h", "--help"))
    {
      std::cerr << help_text << std::endl;
      return 0;
    }
  const double a = cl.follow (double (0.0), "-a");
  const double b = cl.follow (double (1.0), "-b");
  const unsigned int nnodes = cl.follow (100, 2, "-n", "--nnodes");
  const unsigned int nel = nnodes - 1;
  const std::string diffusion = cl.follow ("1.0", 2, "-d", "--diffusion");
  const std::string forcing = cl.follow ("1.0", 2, "-f", "--forcing");

  const double tol = 1e-6;
  const unsigned int maxit = 43;
  const unsigned int overlap = 3;
  const double L = b - a;
  const int mpi_size = 3;

  const double L_loc = L / double (mpi_size);
  const double h = L_loc / ceil (double(nel) / double(mpi_size));

  std::vector<double> a_loc (mpi_size);
  std::vector<double> b_loc (mpi_size);
  std::vector<unsigned int> nel_loc (mpi_size);
  std::vector<unsigned int> ndof_loc (mpi_size);
  std::vector<fem_1d<double>*> subproblems (mpi_size);
  
  coeff<double> a_coeff (diffusion);
  coeff<double> f_coeff (forcing);

  int mpi_rank;
  for (mpi_rank = 0; mpi_rank < mpi_size; ++mpi_rank)
    {
      a_loc[mpi_rank] = a + mpi_rank * L_loc;
      b_loc[mpi_rank] = a_loc[mpi_rank] + L_loc;
      nel_loc[mpi_rank] = ceil (double(nel) / double(mpi_size));
      if (mpi_rank > 0)
        {
          a_loc[mpi_rank] -= overlap * h;
          nel_loc[mpi_rank] += overlap;
        }
      if (mpi_rank < mpi_size - 1)
        {
          b_loc[mpi_rank] += overlap * h;
          nel_loc[mpi_rank] += overlap;
        }
      ndof_loc[mpi_rank] = nel_loc[mpi_rank] + 1;
      fem_1d<double>* tmp = new fem_1d<double>(new mesh<double>
                                               (a_loc[mpi_rank],
                                                b_loc[mpi_rank],
                                                ndof_loc[mpi_rank]));
      subproblems[mpi_rank] = tmp;
      
      subproblems[mpi_rank]->set_diffusion_coefficient (a_coeff);
      subproblems[mpi_rank]->set_source_coefficient (f_coeff);
      subproblems[mpi_rank]->assemble ();
      subproblems[mpi_rank]->set_dirichlet (fem_1d<double>::left_boundary, 0.0);
      subproblems[mpi_rank]->set_dirichlet (fem_1d<double>::right_boundary, 0.0);
      subproblems[mpi_rank]->solve ();
    };

  for (unsigned int it = 0; it < maxit; ++it)
    for (mpi_rank = 0; mpi_rank < mpi_size; ++mpi_rank)
      {
        if (mpi_rank > 0)
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::left_boundary,
             subproblems[mpi_rank-1]->result ()
             [ndof_loc[mpi_rank-1]-1-2*overlap]);
        else
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::left_boundary, 0.0);
        if (mpi_rank < mpi_size - 1)
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::right_boundary,
             subproblems[mpi_rank+1]->result ()
             [2*overlap]);
        else
          subproblems[mpi_rank]->set_dirichlet
            (fem_1d<double>::right_boundary, 0.0);
        subproblems[mpi_rank]->solve ();
      }

  for (mpi_rank = 0; mpi_rank < mpi_size; ++mpi_rank)
    for (unsigned int ii = 0; ii < ndof_loc[mpi_rank]; ++ii)
      std::cout << subproblems[mpi_rank]->m->nodes[ii] << " "
                << subproblems[mpi_rank]->result ()(ii, 0)
                << std::endl;

  return 0;
};

Int
main ( Int argc, char** argv )
{
    //Setup main communicator
    boost::shared_ptr< Epetra_Comm > comm;

    //Setup MPI variables
    Int numberOfProcesses (1);
    Int rank (0);

#ifdef HAVE_MPI
    MPI_Init ( &argc, &argv );

    MPI_Comm_size ( MPI_COMM_WORLD, &numberOfProcesses );
    MPI_Comm_rank ( MPI_COMM_WORLD, &rank );
#endif

    if ( rank == 0 )
    {
        std::cout << std::endl;
        std::cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << std::endl;
        std::cout << " THE ZERO DIMENSIONAL SOLVER IS AN ALPHA VERSION UNDER STRONG DEVELOPMENT" << std::endl;
        std::cout << "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$" << std::endl << std::endl;

        std::cout << "MPI Processes: " << numberOfProcesses << std::endl;
    }

#ifdef HAVE_MPI
    if ( numberOfProcesses > 1 )
    {
        if ( rank == 0 )
        {
            std::cout << "test_ZeroDimensional not enabled in parallel, failing gracefully." << std::endl;
            std::cout << "MPI Finalization" << std::endl;
        }
        MPI_Finalize();
        return EXIT_FAILURE;
    }
#endif

#ifdef EPETRA_MPI
    if ( rank == 0 )
    {
        std::cout << "MPI Epetra Initialization ... " << std::endl;
    }
    comm.reset ( new Epetra_MpiComm ( MPI_COMM_WORLD ) );
#else
    std::cout << "SERIAL Epetra Initialization ... " << std::endl;
    comm.reset ( new Epetra_SerialComm() );
#endif

    bool exitFlag = EXIT_SUCCESS;

#if ( defined(HAVE_NOX_THYRA) && defined(HAVE_TRILINOS_RYTHMOS) )
    // Command line parameters
    GetPot commandLine ( argc, argv );
    const bool check = commandLine.search ( 2, "-c", "--check" );
    string fileName  = commandLine.follow ( "data", 2, "-f", "--file" );

    // SetupData
    GetPot dataFile ( fileName  + ".dat" );

    std::string circuitDataFile = dataFile ( "0D_Model/CircuitDataFile", "./inputFile.dat" );
    BCInterface0D< ZeroDimensionalBCHandler, ZeroDimensionalData >  zeroDimensionalBC;
    zeroDimensionalBC.createHandler();
    zeroDimensionalBC.fillHandler ( circuitDataFile, "Files" );

    boost::shared_ptr< ZeroDimensionalData > zeroDimensionalData ( new ZeroDimensionalData );
    zeroDimensionalData->setup ( dataFile, zeroDimensionalBC.handler() );

    boost::shared_ptr< ZeroDimensionalSolver > zeroDimensionalSolver ( new ZeroDimensionalSolver ( zeroDimensionalData->unknownCounter(), comm, zeroDimensionalData->circuitData() ) );
    zeroDimensionalSolver->setup ( zeroDimensionalData->solverData() );

    zeroDimensionalData->showMe();

    // SetupModel
    zeroDimensionalData->dataTime()->setInitialTime (0);
    zeroDimensionalData->initializeSolution();

    // Create output folder
    if ( comm->MyPID() == 0 )
    {
        mkdir ( "output", 0777 );
    }

    // Save initial solution
    zeroDimensionalData->saveSolution();

    zeroDimensionalData->dataTime()->updateTime();
    zeroDimensionalData->dataTime()->setInitialTime (zeroDimensionalData->dataTime()->time() );

    // Definitions for the time loop
    LifeChrono chronoTotal;
    LifeChrono chronoSystem;
    LifeChrono chronoIteration;

    Int count = 0;
    chronoTotal.start();

    for ( ; zeroDimensionalData->dataTime()->canAdvance() ; zeroDimensionalData->dataTime()->updateTime(), ++count )
//.........这里部分代码省略.........