C++ NumericMatrix::cols方法代码示例

// LMM scan of a single chromosome with interactive covariates
// this version should be fast but requires more memory
// (since we first expand the genotype probabilities to probs x intcovar)
// and this one allows weights for the individuals (the same for all phenotypes)
//
// genoprobs = 3d array of genotype probabilities (individuals x genotypes x positions)
// pheno     = matrix with one column of numeric phenotypes
//             (no missing data allowed)
// addcovar  = additive covariates (an intercept, at least)
// intcovar  = interactive covariates (should also be included in addcovar)
// eigenvec  = matrix of transposed eigenvectors of variance matrix
// weights   = vector of weights (really the SQUARE ROOT of the weights)
//
// output    = vector of log likelihood values
//
// [[Rcpp::export]]
NumericVector scan_pg_onechr_intcovar_highmem(const NumericVector& genoprobs,
                                              const NumericMatrix& pheno,
                                              const NumericMatrix& addcovar,
                                              const NumericMatrix& intcovar,
                                              const NumericMatrix& eigenvec,
                                              const NumericVector& weights,
                                              const double tol=1e-12)
{
    const unsigned int n_ind = pheno.rows();
    if(pheno.cols() != 1)
        throw std::range_error("ncol(pheno) != 1");
    const Dimension d = genoprobs.attr("dim");
    const unsigned int n_pos = d[2];
    if(n_ind != d[0])
        throw std::range_error("nrow(pheno) != nrow(genoprobs)");
    if(n_ind != addcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(addcovar)");
    if(n_ind != intcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(intcovar)");
    if(n_ind != weights.size())
        throw std::range_error("nrow(pheno) != length(weights)");
    if(n_ind != eigenvec.rows())
        throw std::range_error("ncol(pheno) != nrow(eigenvec)");
    if(n_ind != eigenvec.cols())
        throw std::range_error("ncol(pheno) != ncol(eigenvec)");

    // expand genotype probabilities to include geno x interactive covariate
    NumericVector genoprobs_rev = expand_genoprobs_intcovar(genoprobs, intcovar);

    // pre-multiply everything by eigenvectors
    genoprobs_rev = matrix_x_3darray(eigenvec, genoprobs_rev);
    NumericMatrix addcovar_rev = matrix_x_matrix(eigenvec, addcovar);
    NumericMatrix pheno_rev = matrix_x_matrix(eigenvec, pheno);

    // multiply everything by the (square root) of the weights
    // (weights should ALREADY be the square-root of the real weights)
    addcovar_rev = weighted_matrix(addcovar_rev, weights);
    pheno_rev = weighted_matrix(pheno_rev, weights);
    genoprobs_rev = weighted_3darray(genoprobs_rev, weights);

    // regress out the additive covariates
    genoprobs_rev = calc_resid_linreg_3d(addcovar_rev, genoprobs_rev, tol);
    pheno_rev = calc_resid_linreg(addcovar_rev, pheno_rev, tol);

    // now the scan, return RSS
    NumericMatrix rss = scan_hk_onechr_nocovar(genoprobs_rev, pheno_rev, tol);

    // 0.5*sum(log(weights)) [since these are sqrt(weights)]
    double sum_logweights = sum(log(weights));

    // calculate log likelihood
    NumericVector result(n_pos);
    for(unsigned int pos=0; pos<n_pos; pos++)
        result[pos] = -(double)n_ind/2.0*log(rss[pos]) + sum_logweights;

    return result;
}

// LMM scan of a single chromosome with interactive covariates
// this version uses less memory but will be slower
// (since we need to work with each position, one at a time)
// and this one allows weights for the individuals (the same for all phenotypes)
//
// genoprobs = 3d array of genotype probabilities (individuals x genotypes x positions)
// pheno     = matrix with one column of numeric phenotypes
//             (no missing data allowed)
// addcovar  = additive covariates (an intercept, at least)
// intcovar  = interactive covariates (should also be included in addcovar)
// eigenvec  = matrix of transposed eigenvectors of variance matrix
// weights   = vector of weights (really the SQUARE ROOT of the weights)
//
// output    = vector of log likelihood values
//
// [[Rcpp::export]]
NumericVector scan_pg_onechr_intcovar_lowmem(const NumericVector& genoprobs,
                                             const NumericMatrix& pheno,
                                             const NumericMatrix& addcovar,
                                             const NumericMatrix& intcovar,
                                             const NumericMatrix& eigenvec,
                                             const NumericVector& weights,
                                             const double tol=1e-12)
{
    const unsigned int n_ind = pheno.rows();
    if(pheno.cols() != 1)
        throw std::range_error("ncol(pheno) != 1");
    const Dimension d = genoprobs.attr("dim");
    const unsigned int n_pos = d[2];
    if(n_ind != d[0])
        throw std::range_error("nrow(pheno) != nrow(genoprobs)");
    if(n_ind != addcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(addcovar)");
    if(n_ind != intcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(intcovar)");
    if(n_ind != weights.size())
        throw std::range_error("ncol(pheno) != length(weights)");
    if(n_ind != eigenvec.rows())
        throw std::range_error("ncol(pheno) != nrow(eigenvec)");
    if(n_ind != eigenvec.cols())
        throw std::range_error("ncol(pheno) != ncol(eigenvec)");

    NumericVector result(n_pos);

    // multiply pheno by eigenvectors and then by weights
    NumericMatrix pheno_rev = matrix_x_matrix(eigenvec, pheno);
    pheno_rev = weighted_matrix(pheno_rev, weights);

    // 0.5*sum(log(weights)) [since these are sqrt(weights)]
    double sum_logweights = sum(log(weights));

    for(unsigned int pos=0; pos<n_pos; pos++) {
        Rcpp::checkUserInterrupt();  // check for ^C from user

        // form X matrix
        NumericMatrix X = formX_intcovar(genoprobs, addcovar, intcovar, pos, true);

        // multiply by eigenvectors
        X = matrix_x_matrix(eigenvec, X);

        // multiply by weights
        X = weighted_matrix(X, weights);

        // do regression
        NumericVector rss = calc_rss_linreg(X, pheno_rev, tol);
        result[pos] = -(double)n_ind/2.0*log(rss[0]) + sum_logweights;
    }

    return result;
}

// [[Rcpp::export]]
SEXP HandMadeKalmanFilterConstantCoeffCpp(NumericVector a0
  , NumericMatrix P0, NumericMatrix T, NumericMatrix Z , NumericMatrix HH, NumericMatrix GG, NumericMatrix yt)
{
  // convert data to Eigen class
  Eigen::Map<Eigen::VectorXd > a0e(&a0[0], (size_t)(a0.size()));
  Eigen::Map<Eigen::VectorXd > P0e(&P0[0], P0.rows(), P0.cols());
  Eigen::Map<Eigen::MatrixXd > Te(&T[0], T.rows(), T.cols());
  Eigen::Map<Eigen::MatrixXd > Ze(&Z[0], Z.rows(), Z.cols());
  Eigen::Map<Eigen::MatrixXd > HHe(&HH[0], HH.rows(), HH.cols());
  Eigen::Map<Eigen::MatrixXd > GGe(&GG[0], GG.rows(), GG.cols());
  
  // HMKF initialization block
  cout << T.rows() << " " << Z.rows() << " " << HH.rows() << " " << GG.rows() << endl;
  HMKF kf = HMKF(T.rows(), Z.rows());
  kf.setModel(Te, Ze, HHe, GGe);
  kf.initialize(a0e, P0e);

  // filtering steps
  NumericVector x(yt.cols()), xp(yt.cols()), V(yt.cols());
  for(int i=0; i!=yt.cols(); ++i){
    kf.predict();
    NumericMatrix::Column col = yt(_,i);
//    std::cout << "y:" << col[0] << std::endl;
    kf.filter(Map<Eigen::VectorXd >(&col[0], col.size()));
//    xp[i]= (kf.getxp())(0);
    x[i] = (kf.getx())(0);
//    V[i] = (kf.getV())(0,0);
  }
  return List::create(_("x") = x, _("xp") = xp, _("V") = V);
}

// Scan a single chromosome with additive covariates and weights
//
// genoprobs = 3d array of genotype probabilities (individuals x genotypes x positions)
// pheno     = matrix of numeric phenotypes (individuals x phenotypes)
//             (no missing data allowed, values should be in [0,1])
// addcovar  = additive covariates (an intercept, at least)
// weights   = vector of weights
//
// output    = matrix of (weighted) residual sums of squares (RSS) (phenotypes x positions)
//
// [[Rcpp::export]]
NumericMatrix scan_binary_onechr_weighted(const NumericVector& genoprobs,
                                          const NumericMatrix& pheno,
                                          const NumericMatrix& addcovar,
                                          const NumericVector& weights,
                                          const int maxit=100,
                                          const double tol=1e-6,
                                          const double qr_tol=1e-12,
                                          const double eta_max=30.0)
{
    const int n_ind = pheno.rows();
    if(Rf_isNull(genoprobs.attr("dim")))
        throw std::invalid_argument("genoprobs should be a 3d array but has no dim attribute");
    const Dimension d = genoprobs.attr("dim");
    if(d.size() != 3)
        throw std::invalid_argument("genoprobs should be a 3d array");
    if(n_ind != d[0])
        throw std::range_error("nrow(pheno) != nrow(genoprobs)");
    if(n_ind != addcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(addcovar)");
    if(n_ind != weights.size())
        throw std::range_error("nrow(pheno) != length(weights)");
    const int n_pos = d[2];
    const int n_gen = d[1];
    const int n_add = addcovar.cols();
    const int g_size = n_ind * n_gen;
    const int n_phe = pheno.cols();

    NumericMatrix result(n_phe, n_pos);
    NumericMatrix X(n_ind, n_gen+n_add);

    if(n_add > 0) // paste in covariates, if present
        std::copy(addcovar.begin(), addcovar.end(), X.begin() + g_size);

    for(int pos=0, offset=0; pos<n_pos; pos++, offset += g_size) {
        Rcpp::checkUserInterrupt();  // check for ^C from user

        // copy genoprobs for this pos into a matrix
        std::copy(genoprobs.begin() + offset, genoprobs.begin() + offset + g_size, X.begin());

        for(int phe=0; phe<n_phe; phe++) {
            // calc rss and paste into ith column of result
            result(phe,pos) = calc_ll_binreg_weighted(X, pheno(_,phe), weights, maxit, tol, qr_tol, eta_max);
        }
    }

    return result;
}

// calculate distance between two lists of points
//
// pt1 and pt2 are each matrices with two columns
// [[Rcpp::export]]
NumericMatrix calc_dist_pt2pt(NumericMatrix pt1, NumericMatrix pt2)
{
    int nr1 = pt1.rows(), nr2 = pt2.rows();
    int nc1 = pt1.cols(), nc2 = pt2.cols();
    if(nc1 != 2) throw std::invalid_argument("pt1 should have two columns");
    if(nc2 != 2) throw std::invalid_argument("pt2 should have two columns");

    NumericMatrix result(nr1,nr2);

    for(int i=0; i<nr1; i++) {
        for(int j=0; j<nr2; j++) {
            result(i,j) = calc_dist_pt2pt_one(pt1(i,_), pt2(j,_));
        }
    }

    return result;
}

// [[Rcpp::export]]
NumericMatrix term_score(NumericMatrix beta) {
  // calculate term score (Blei and Lafferty 2009)
  int W = beta.rows();
  int K = beta.cols();
  NumericMatrix term_score(W, K);
  
  for (int w = 0; w < W; w++) {
    double product = 0.0;
    for (int k = 0; k < K; k++) {
      product *= beta(w, k);
    }
    for (int k = 0; k < K; k++) {
      term_score(w, k) = beta(w, k) * log(beta(w, k) / pow(product, (1 / K)));
    }
  }
  
  return term_score;
}

// [[Rcpp::export]]
NumericVector row_medians(NumericMatrix toSort) {
  int n = toSort.rows();
  int medN  = toSort.cols();
  NumericVector meds = NumericVector(n);
  for (int i = 0; i < n; i++) {
    NumericVector curRow = toSort.row(i);
    std::nth_element(curRow.begin(), curRow.begin() + curRow.size()/2 - 1, curRow.end());
    double med1 = curRow[curRow.size()/2 - 1];
    if (medN % 2 == 0) {
      std::nth_element(curRow.begin(), curRow.begin() + curRow.size()/2, curRow.end());
      double med2 = curRow[curRow.size()/2];
      meds[i] = (med1 + med2)/2.0;
    } else {
      meds[i] = med1;
    }
  }

  return meds;
}

// Scan a single chromosome with interactive covariates
// this version uses less memory but will be slower
// (since we need to work with each position, one at a time)
// and this one allows weights for the individuals (the same for all phenotypes)
//
// genoprobs = 3d array of genotype probabilities (individuals x genotypes x positions)
// pheno     = matrix of numeric phenotypes (individuals x phenotypes)
//             (no missing data allowed)
// addcovar  = additive covariates (an intercept, at least)
// intcovar  = interactive covariates (should also be included in addcovar)
// weights   = vector of weights
//
// output    = matrix of residual sums of squares (RSS) (phenotypes x positions)
//
// [[Rcpp::export]]
NumericMatrix scan_binary_onechr_intcovar_weighted_lowmem(const NumericVector& genoprobs,
                                                          const NumericMatrix& pheno,
                                                          const NumericMatrix& addcovar,
                                                          const NumericMatrix& intcovar,
                                                          const NumericVector& weights,
                                                          const int maxit=100,
                                                          const double tol=1e-6,
                                                          const double qr_tol=1e-12,
                                                          const double eta_max=30.0)
{
    const int n_ind = pheno.rows();
    if(Rf_isNull(genoprobs.attr("dim")))
        throw std::invalid_argument("genoprobs should be a 3d array but has no dim attribute");
    const Dimension d = genoprobs.attr("dim");
    if(d.size() != 3)
        throw std::invalid_argument("genoprobs should be a 3d array");
    const int n_pos = d[2];
    const int n_phe = pheno.cols();
    if(n_ind != d[0])
        throw std::range_error("nrow(pheno) != nrow(genoprobs)");
    if(n_ind != addcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(addcovar)");
    if(n_ind != intcovar.rows())
        throw std::range_error("nrow(pheno) != nrow(intcovar)");

    NumericMatrix result(n_phe, n_pos);

    for(int pos=0; pos<n_pos; pos++) {
        Rcpp::checkUserInterrupt();  // check for ^C from user

        // form X matrix
        NumericMatrix X = formX_intcovar(genoprobs, addcovar, intcovar, pos, true);

        for(int phe=0; phe<n_phe; phe++) {
            // do regression
            result(phe,pos) = calc_ll_binreg_weighted(X, pheno(_,phe), weights, maxit, tol, qr_tol, eta_max);
        }
    }

    return result;
}

// [[Rcpp::export]]
NumericVector calc_coherence(List ctot, NumericMatrix top_words_topics) {
  int K = top_words_topics.cols();
  int M = top_words_topics.rows();
  NumericMatrix cond_topic_count = ctot["cond_topic_count"];
  NumericVector coherence_scores(K);
  NumericVector composer_id = ctot["composer_id"];
  for (int k = 0; k < K; k++) {
    for (int m = 2; m < M; m++) {
      for (int l = 0; l < m-1; l++) {
        double doc_freq = 0.0;
        double co_doc_freq = 0.0;
        for (int i = 0; i < composer_id.size(); i++) {
          doc_freq += (double) (composer_id[i] == top_words_topics(l, k));
          co_doc_freq += (double) (composer_id[i] == top_words_topics(m, k) &&
                                   composer_id[i] == top_words_topics(l, k));
        }
        coherence_scores[k] += log((co_doc_freq + 1.0) / doc_freq);
      }
    }
  }
  return coherence_scores;
}

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();
//.........这里部分代码省略.........

RcppExport SEXP robustMatrixTransform( SEXP r_matrix )
{
try
{
  typedef double RealType;
  NumericMatrix M = as< NumericMatrix >( r_matrix );
  NumericMatrix outMat( M.rows(), M.cols() );
  unsigned long rows = M.rows();
  for( unsigned long j = 0; j < M.cols(); j++ )
    {
    NumericVector Mvec = M(_, j);
    NumericVector rank = M(_, j);
    for( unsigned int i = 0; i < rows; i++ )
      {
      RealType   rankval = 0;
      RealType   xi = Mvec(i);
      for( unsigned int k = 0; k < rows; k++ )
        {
        RealType yi = Mvec(k);
        RealType diff = fabs(xi - yi);
        if( diff > 0 )
          {
          RealType val = (xi - yi) / diff;
          rankval += val;
          }
        }
      rank(i) = rankval / rows;
      }
    outMat(_, j) = rank;
    }
  // this passes rcpp data to vnl matrix ... safely?
  if ( 1 == 0 )
    {
    // Further notes on this:
    // 1. see multiChannel cpp for how to pass image list
    // 2. separate implementation for eanat and sccan
    // 3. deal with output as only matrix?
    // 4. .... ?
    std::vector<RealType> dat =
      Rcpp::as< std::vector<RealType> >( outMat );
    const double* _data = &dat[0];
    vnl_matrix<RealType> vnlmat( _data , M.cols(), M.rows()  );
    vnlmat = vnlmat.transpose();
    }
  return wrap( outMat );
}
catch( itk::ExceptionObject & err )
  {
  Rcpp::Rcout << "ITK ExceptionObject caught !" << std::endl;
  forward_exception_to_r( err );
  }
catch( const std::exception& exc )
  {
  Rcpp::Rcout << "STD ExceptionObject caught !" << std::endl;
  forward_exception_to_r( exc );
  }
catch(...)
  {
	Rcpp::stop("c++ exception (unknown reason)");
  }
return Rcpp::wrap(NA_REAL); //not reached
}