C++ Pointer::GetPixel方法代码示例

SEXP addNeighborhoodToImageHelper(
  SEXP r_antsimage,
  SEXP r_center,
  SEXP r_rad,
  SEXP r_vec)
{
  typedef typename ImageType::Pointer  ImagePointerType;
  const unsigned int ImageDimension = ImageType::ImageDimension;
  typedef float                        PixelType;
  typename ImageType::Pointer image =
    Rcpp::as< ImagePointerType >( r_antsimage );
  Rcpp::NumericVector center( r_center );
  Rcpp::NumericVector rad( r_rad );
  Rcpp::NumericVector intvec( r_vec );
  if ( center.size() != ImageDimension )
    Rcpp::stop("addNeighborhoodToImageHelper dim error.");
  typename itk::NeighborhoodIterator<ImageType>::SizeType nSize;
  typename ImageType::IndexType ind;
  ind.Fill( 0 );
  for ( unsigned int i=0; i<ImageDimension; i++ )
    {
    nSize[i] = rad[i];
    ind[i] = center[i]; // R coords to ITK
    }
  itk::NeighborhoodIterator<ImageType> nit( nSize, image,
    image->GetLargestPossibleRegion() ) ;
// for each location in nitSearch, compute the correlation
// of the intvec with the nit neighborhood
  nit.SetLocation( ind );
  for( unsigned int i = 0; i < nit.Size(); i++ )
    {
    typename ImageType::IndexType ind2 = nit.GetIndex(i);
    PixelType lval = image->GetPixel( ind2 );
    image->SetPixel( ind2, lval + intvec[i] );
    }
  return 0;
}

SEXP eigenanatomyCppHelper(
  NumericMatrix X,
  SEXP r_mask,
  RealType sparseness,
  IntType nvecs,
  IntType its,
  IntType cthresh,
  RealType z,
  RealType smooth,
//  NumericMatrix initializationMatrix,
  Rcpp::List initializationList,
  IntType covering,
  RealType ell1,
  IntType verbose,
  IntType powerit,
  RealType priorWeight )
{
  enum { Dimension = ImageType::ImageDimension };
  typename ImageType::RegionType region;
  typedef typename ImageType::PixelType PixelType;
  typedef typename ImageType::Pointer ImagePointerType;
  typedef double                                        Scalar;
  typedef itk::ants::antsSCCANObject<ImageType, Scalar> SCCANType;
  typedef typename SCCANType::MatrixType                vMatrix;
  typename SCCANType::Pointer sccanobj = SCCANType::New();

  typename ImageType::Pointer mask = Rcpp::as<ImagePointerType>( r_mask );
  bool maskisnull = mask.IsNull();
// deal with the initializationList, if any
  unsigned int nImages = initializationList.size();
  if ( ( nImages > 0 ) && ( !maskisnull ) )
    {
    itk::ImageRegionIteratorWithIndex<ImageType> it( mask,
      mask->GetLargestPossibleRegion() );
    vMatrix priorROIMat( nImages , X.cols() );
    priorROIMat.fill( 0 );
    for ( unsigned int i = 0; i < nImages; i++ )
      {
      typename ImageType::Pointer init =
        Rcpp::as<ImagePointerType>( initializationList[i] );
      unsigned long ct = 0;
      it.GoToBegin();
      while ( !it.IsAtEnd() )
        {
        PixelType pix = it.Get();
        if ( pix >= 0.5 )
          {
          pix = init->GetPixel( it.GetIndex() );
          priorROIMat( i, ct ) = pix;
          ct++;
          }
        ++it;
        }
      }
    sccanobj->SetMatrixPriorROI( priorROIMat );
    nvecs = nImages;
    }
  sccanobj->SetPriorWeight( priorWeight );
  sccanobj->SetLambda( priorWeight );
// cast hack from Rcpp type to sccan type
  std::vector<double> xdat =
      Rcpp::as< std::vector<double> >( X );
  const double* _xdata = &xdat[0];
  vMatrix vnlX( _xdata , X.cols(), X.rows()  );
  vnlX = vnlX.transpose();

  sccanobj->SetGetSmall( false  );
  vMatrix priorROIMat;

//    sccanobj->SetMatrixPriorROI( priorROIMat);
//    sccanobj->SetMatrixPriorROI2( priorROIMat );
  sccanobj->SetCovering( covering );
  sccanobj->SetSilent(  ! verbose  );
  if( ell1 > 0 )
    {
    sccanobj->SetUseL1( true );
    }
  else
    {
    sccanobj->SetUseL1( false );
    }
  sccanobj->SetGradStep( vnl_math_abs( ell1 ) );
  sccanobj->SetMaximumNumberOfIterations( its );
  sccanobj->SetRowSparseness( z );
  sccanobj->SetSmoother( smooth );
  if ( sparseness < 0 ) sccanobj->SetKeepPositiveP(false);
  sccanobj->SetSCCANFormulation(  SCCANType::PQ );
  sccanobj->SetFractionNonZeroP( fabs( sparseness ) );
  sccanobj->SetMinClusterSizeP( cthresh );
  sccanobj->SetMatrixP( vnlX );
//  sccanobj->SetMatrixR( r ); // FIXME
  sccanobj->SetMaskImageP( mask );
  RealType truecorr = 0;
  if( powerit == 1 )
    {
    truecorr = sccanobj->SparseReconHome( nvecs );
    }
  else if ( priorWeight > 1.e-12 )
    truecorr = sccanobj->SparseReconPrior( nvecs, true );
  else truecorr = sccanobj->SparseRecon(nvecs);
//.........这里部分代码省略.........

SEXP sccanCppHelper(
  NumericMatrix X,
  NumericMatrix Y,
  SEXP r_maskx,
  SEXP r_masky,
  RealType sparsenessx,
  RealType sparsenessy,
  IntType nvecs,
  IntType its,
  IntType cthreshx,
  IntType cthreshy,
  RealType z,
  RealType smooth,
  Rcpp::List initializationListx,
  Rcpp::List initializationListy,
  IntType covering,
  RealType ell1,
  IntType verbose,
  RealType priorWeight )
{
  enum { Dimension = ImageType::ImageDimension };
  typename ImageType::RegionType region;
  typedef typename ImageType::PixelType PixelType;
  typedef typename ImageType::Pointer ImagePointerType;
  typedef double                                        Scalar;
  typedef itk::ants::antsSCCANObject<ImageType, Scalar> SCCANType;
  typedef typename SCCANType::MatrixType                vMatrix;
  typename SCCANType::Pointer sccanobj = SCCANType::New();

  typename ImageType::Pointer maskx = Rcpp::as<ImagePointerType>( r_maskx );
  typename ImageType::Pointer masky = Rcpp::as<ImagePointerType>( r_masky );

  bool maskxisnull = maskx.IsNull();
  bool maskyisnull = masky.IsNull();
// deal with the initializationList, if any
  unsigned int nImagesx = initializationListx.size();
  if ( ( nImagesx > 0 ) && ( !maskxisnull ) )
    {
    itk::ImageRegionIteratorWithIndex<ImageType> it( maskx,
      maskx->GetLargestPossibleRegion() );
    vMatrix priorROIMatx( nImagesx , X.cols() );
    priorROIMatx.fill( 0 );
    for ( unsigned int i = 0; i < nImagesx; i++ )
      {
      typename ImageType::Pointer init =
        Rcpp::as<ImagePointerType>( initializationListx[i] );
      unsigned long ct = 0;
      it.GoToBegin();
      while ( !it.IsAtEnd() )
        {
        PixelType pix = it.Get();
        if ( pix >= 0.5 )
          {
          pix = init->GetPixel( it.GetIndex() );
          priorROIMatx( i, ct ) = pix;
          ct++;
          }
        ++it;
        }
      }
    sccanobj->SetMatrixPriorROI( priorROIMatx );
    nvecs = nImagesx;
    }
  unsigned int nImagesy = initializationListy.size();
  if ( ( nImagesy > 0 ) && ( !maskyisnull ) )
    {
    itk::ImageRegionIteratorWithIndex<ImageType> it( masky,
      masky->GetLargestPossibleRegion() );
    vMatrix priorROIMaty( nImagesy , Y.cols() );
    priorROIMaty.fill( 0 );
    for ( unsigned int i = 0; i < nImagesy; i++ )
      {
      typename ImageType::Pointer init =
        Rcpp::as<ImagePointerType>( initializationListy[i] );
      unsigned long ct = 0;
      it.GoToBegin();
      while ( !it.IsAtEnd() )
        {
        PixelType pix = it.Get();
        if ( pix >= 0.5 )
          {
          pix = init->GetPixel( it.GetIndex() );
          priorROIMaty( i, ct ) = pix;
          ct++;
          }
        ++it;
        }
      }
    sccanobj->SetMatrixPriorROI2( priorROIMaty );
    nvecs = nImagesy;
    }

  sccanobj->SetPriorWeight( priorWeight );
  sccanobj->SetLambda( priorWeight );
// cast hack from Rcpp type to sccan type
  std::vector<double> xdat =
      Rcpp::as< std::vector<double> >( X );
  const double* _xdata = &xdat[0];
  vMatrix vnlX( _xdata , X.cols(), X.rows()  );
  vnlX = vnlX.transpose();
//.........这里部分代码省略.........

SEXP jointLabelFusionNeighborhoodSearchHelper(
  SEXP r_intvec,
  SEXP r_center,
  unsigned int rad,
  unsigned int radSearch,
  SEXP r_antsimage,
  SEXP r_antsimageseg)
{
  unsigned int segval = 0;
  typedef typename ImageType::Pointer  ImagePointerType;
  const unsigned int ImageDimension = ImageType::ImageDimension;
  typedef float                        PixelType;
  typename ImageType::Pointer image =
    Rcpp::as< ImagePointerType >( r_antsimage );
  typename ImageType::Pointer imageseg =
      Rcpp::as< ImagePointerType >( r_antsimageseg );
  Rcpp::NumericVector intvec( r_intvec );
  Rcpp::NumericVector outvec =
    Rcpp::NumericVector( intvec.size(), Rcpp::NumericVector::get_na() );
  Rcpp::NumericVector bestvec =
    Rcpp::NumericVector( intvec.size(), Rcpp::NumericVector::get_na() );
  Rcpp::NumericVector outsegvec =
    Rcpp::NumericVector( intvec.size(), Rcpp::NumericVector::get_na() );
  Rcpp::NumericVector bestsegvec =
    Rcpp::NumericVector( intvec.size(), Rcpp::NumericVector::get_na() );
  Rcpp::NumericVector center( r_center );
  if ( center.size() != ImageDimension )
    Rcpp::stop("jointLabelFusionNeighborhoodSearchHelper dim error.");
  typename itk::NeighborhoodIterator<ImageType>::SizeType nSize;
  typename itk::NeighborhoodIterator<ImageType>::SizeType nSizeSearch;
  typename ImageType::IndexType ind;
  ind.Fill( 0 );
  for ( unsigned int i=0; i<ImageDimension; i++ )
    {
    nSize[i] = rad;
    nSizeSearch[i] = radSearch;
    ind[i] = center[i]; // R coords to ITK
    }
  itk::NeighborhoodIterator<ImageType> nit( nSize, image,
    image->GetLargestPossibleRegion() ) ;
  itk::NeighborhoodIterator<ImageType> nitSearch( nSizeSearch, image,
    image->GetLargestPossibleRegion() ) ;
// for each location in nitSearch, compute the correlation
// of the intvec with the nit neighborhood
  nitSearch.SetLocation( ind );
  PixelType bestcor = 1.e11;
  PixelType bestsd = 0;
  PixelType bestmean = 0;
  for( unsigned int i = 0; i < nitSearch.Size(); i++ )
    {
    typename ImageType::IndexType ind2 = nitSearch.GetIndex(i);
    nit.SetLocation( ind2 );
    PixelType outmean = 0;
    PixelType outsd = 0;
    PixelType inmean = 0;
    PixelType insd = 0;
    for ( unsigned int i=0; i < intvec.size(); i++ ) {
      PixelType pix = image->GetPixel( nit.GetIndex(i) );
      outvec[i] = pix;
      outsegvec[i] = imageseg->GetPixel( nit.GetIndex(i) );
      outmean += pix;
      inmean += intvec[i];
      }
    outmean /= ( static_cast<PixelType>(intvec.size()) );
    inmean /= ( static_cast<PixelType>(intvec.size()) );
    for ( unsigned int i=0; i < intvec.size(); i++ ) {
      // should use recursive formula in above loop
      outsd += ( outvec[i] - outmean ) * ( outvec[i] - outmean );
      insd += ( intvec[i] - inmean ) * ( intvec[i] - inmean );
      }
    outsd = sqrt(  outsd  );
    insd = sqrt(  insd  );
    PixelType sum_uv = 0;
    PixelType sum_psearch = 0;
    PixelType ssq_psearch = 0;
    unsigned int n = intvec.size();
    for(unsigned int i = 0; i < n; i++)
      {
      PixelType v = intvec[i];
      PixelType u = outvec[i];
      sum_psearch += u;
      ssq_psearch += u * u;
      sum_uv += u * v;
      }
    PixelType var_u_unnorm = ssq_psearch - sum_psearch * sum_psearch / n;
    if(var_u_unnorm < 1.0e-6)
      var_u_unnorm = 1.0e-6;
    PixelType locor = 0;
    if ( sum_uv > 0 )
      locor = ( -1.0 * (sum_uv * sum_uv) / var_u_unnorm );
    else
      locor = ( sum_uv * sum_uv ) / var_u_unnorm;
//  - (\Sum u_i v_i)^2 / z,   where z = sigma_v^2 * (n-1)
//      locor = locor / ( insd * outsd );
      if ( locor < bestcor )
        {
        segval = imageseg->GetPixel( ind2 );
        for ( unsigned int i=0; i < intvec.size(); i++ ) {
          bestvec[i] = outvec[i];
          bestsegvec[i] = outsegvec[i];
//.........这里部分代码省略.........