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

//f  = f./repmat(sum(f,1),[classN,1]);
vnl_matrix<double> MCLR_SM::Normalize_F_Sum(vnl_matrix<double> f)
{
	// Matrix for normalization
	vnl_matrix<double> norm_matrix(f.rows(),f.cols());
	vnl_vector<double> norm_matrix_row;
	norm_matrix_row.set_size(f.cols());
	
	 //repmat(sum(f,1),[classN,1]);
	for(int i=0;i<f.cols();++i)
	{
      double sum = 0 ;  
	  for(int j=0;j<no_of_classes;++j)
		{
			sum = sum + f(j,i);
		}
		norm_matrix_row(i) = sum;
	}


    for(int i=0;i<no_of_classes;++i)
	{
		norm_matrix.set_row(i,norm_matrix_row); 
	}

// f  = f./repmat(sum(f,1),[classN,1]);
	for(int i=0;i<f.rows();++i)
	{
	  for(int j=0;j<f.cols();++j)
	   {
		  f(i,j) = f(i,j)/norm_matrix(i,j);
	   }
	}

	return f;
}

void pick_indices(const vnl_matrix<double>&dist,
    vnl_matrix<int>&pairs, const double threshold) {
  int m = dist.rows();
  int n = dist.cols();
  vnl_vector<int> row_flag, col_flag;
  col_flag.set_size(n);  col_flag.fill(0);
  row_flag.set_size(n);  row_flag.fill(0);
  std::vector<int> row_index,col_index;
  for (int i = 0; i < m; ++i) {
    double min_dist = dist.get_row(i).min_value();
    if (min_dist < threshold) {
      for (int j = 0; j < n; ++j) {
        if (dist(i,j)==min_dist && col_flag[j] == 0){
          row_index.push_back(i);
          row_flag[i] = 1;
          col_index.push_back(j);
          col_flag[j] = 1;
        }
      }
    }
  }
  pairs.set_size(2, row_index.size());
  for (int i = 0; i<pairs.cols(); ++i){
    pairs(0,i) = row_index[i];
    pairs(1,i) = col_index[i];
  }
}

bool rgrsn_ldp::local_dynamic_programming(const vnl_matrix<double> & probMap, int nNeighborBin,
                                          vcl_vector<int> & optimalBins)
{
    const int N    = probMap.rows();
    const int nBin = probMap.cols();
    
    // dynamic programming
    vnl_matrix<double> accumulatedProbMap = vnl_matrix<double>(N, nBin);
    accumulatedProbMap.fill(0.0);
    vnl_matrix<int> lookbackTable = vnl_matrix<int>(N, nBin);
    lookbackTable.fill(0);
    // copy first row
    for (int c = 0; c<probMap.cols(); c++) {
        accumulatedProbMap[0][c] = probMap[0][c];
        lookbackTable[0][c] = c;
    }
    
    for (int r = 1; r <N; r++) {
        for (int c = 0; c<probMap.cols(); c++) {
            // lookup all possible place in the window
            double max_val = -1;
            int max_index  = -1;
            for (int w = -nNeighborBin; w <= nNeighborBin; w++) {
                if (c + w <0 || c + w >= probMap.cols()) {
                    continue;
                }
                double val = probMap[r][c] + accumulatedProbMap[r-1][c+w];
                if (val > max_val) {
                    max_val = val;
                    max_index = c + w; // most probable path from the [r-1] row, in column c + w
                }
            }
            assert(max_index != -1);
            accumulatedProbMap[r][c] = max_val;
            lookbackTable[r][c]      = max_index;
        }
    }
    
    // lookback the table
    double max_prob    = -1.0;
    int max_prob_index = -1;
    for (int c = 0; c<accumulatedProbMap.cols(); c++) {
        if (accumulatedProbMap[N-1][c] > max_prob) {
            max_prob = accumulatedProbMap[N-1][c];
            max_prob_index = c;
        }
    }
    
    // back track
    optimalBins.push_back(max_prob_index);
    for (int r = N-1; r > 0; r--) {
        int bin = lookbackTable[r][optimalBins.back()];
        optimalBins.push_back(bin);
    }
    assert(optimalBins.size() == N);
    
  //  vcl_reverse(optimalBins.begin(), optimalBins.end());
    return true;
}

bool rgrsn_ldp::dynamic_programming(const vnl_matrix<double> & data, double v_min, double v_max,
                                    unsigned int nBin, int nJumpBin, unsigned int windowSize,
                                    vnl_vector<double> & optimalSignal)
{
    assert(v_min < v_max);
    // raw data to probability map
    const int N = data.rows();
    vnl_matrix<double> probMap = vnl_matrix<double>(N, nBin);
    double interval = (v_max - v_min)/nBin;
    for (int r = 0; r<N; r++) {
        for (int c = 0; c<data.cols(); c++) {
            int num = value_to_bin_number(v_min, interval, data[r][c], nBin);
            probMap[r][num] += 1.0;
        }
    }
    probMap /= data.cols(); // normalize
    
    // save prob
    {
      //  vnl_matlab_filewrite awriter("prob.mat");
      //  awriter.write(probMap, "prob");
     //   printf("save to prob.mat.\n");
    }
  
    vcl_vector<double> optimalValues(N, 0);
    vcl_vector<int> numValues(N, 0);      // multiple values from local dynamic programming
    for (int i = 0; i<=N - windowSize; i++) {
        // get a local probMap;
        vnl_matrix<double> localProbMap = probMap.extract(windowSize, probMap.cols(), i, 0);
        vcl_vector<int> localOptimalBins;
        rgrsn_ldp::local_dynamic_programming(localProbMap, nJumpBin, localOptimalBins);
        assert(localOptimalBins.size() == windowSize);
        for (int j = 0; j < localOptimalBins.size(); j++) {
            numValues[j + i]     += 1;
            optimalValues[j + i] += bin_number_to_value(v_min, interval, localOptimalBins[j]);
        }
        
        // test
        if (0 && i == 0)
        {
            printf("test first window output\n");
            for (int j = 0; j<optimalValues.size() && j<windowSize; j++) {
                printf("%f ", optimalValues[j]);
            }
            printf("\n");
        }
    }
    //
    for (int i = 0; i<optimalValues.size(); i++) {
        optimalValues[i] /= numValues[i];
    }
    optimalSignal = vnl_vector<double>(&optimalValues[0], (int)optimalValues.size());
    return true;
}

/*template <typename T1, typename T2>
itk::Image<T1,2>::Pointer vnlmat_to_itkimage(vnl_matrix<T2> mat,float scale = 1.0)
{
	typedef typename itk::Image<T1,2> ImageType;
	ImageType::SizeType size;
	ImageType::IndexType index;
	ImageType::RegionType region;
	index.Fill(0);
	size[0] = mat.cols();
	size[1] = mat.rows();
	region.SetSize(size);
	region.SetIndex(index);
	ImageType::Pointer im = ImageType::New();
	im->SetRegions(region);
	im->Allocate();
	for(int xco = 0; xco<size[1]; xco++)
	{
		for(int yco = 0; yco < size[0]; yco++)
		{
			index[1] = xco;
			index[0] = yco;
			im->SetPixel(index,mat[xco][yco]*scale);
		}
	}
	return im;

}
*/

vnl_matrix<float> cpx_to_abs(vnl_matrix< std::complex< float> > mat)
{
	vnl_matrix<float> out(mat.rows(),mat.cols());

	for(int x= 0; x< mat.rows(); x++)
	{
		for(int y = 0; y < mat.cols(); y++)
		{
			out[x][y] = abs(mat[x][y]);
		}
	}
	return out;

bool rgrsn_ldp::local_viterbi(const vnl_matrix<double> & data,
                              double resolution,
                              const vnl_vector<double> & transition,
                              unsigned int window_size,
                              vnl_vector<double> & optimal_signal)
{
    assert(resolution > 0.0);
    assert(transition.size()%2 == 1);
    
    const double min_v = data.min_value();
    const double max_v = data.max_value();
    const int nBin = (max_v - min_v)/resolution;
    
    // raw data to probability map
    // quantilization
    const int N = data.rows();
    vnl_matrix<double> probMap = vnl_matrix<double>(N, nBin);
    for (int r = 0; r<N; r++) {
        for (int c = 0; c<data.cols(); c++) {
            int num = value_to_bin_number(min_v, resolution, data[r][c], nBin);
            probMap[r][num] += 1.0;
        }
    }
    probMap /= data.cols(); // normalization
    
    vcl_vector<double> optimalValues(N, 0);
    vcl_vector<int> numValues(N, 0);      // multiple values from local dynamic programming
    
    for (int i = 0; i <= N - window_size; i++) {
        // get a local probMap;
        vnl_matrix<double> localProbMap = probMap.extract(window_size, probMap.cols(), i, 0);
        vcl_vector<int> localOptimalBins;
        rgrsn_ldp::viterbi(localProbMap, transition, localOptimalBins);
        assert(localOptimalBins.size() == window_size);
        for (int j = 0; j < localOptimalBins.size(); j++) {
            double value = bin_number_to_value(min_v, resolution, localOptimalBins[j]);
            numValues[j + i]     += 1;
            optimalValues[j + i] += value;
        }
    }
    
    // average all optimal path as final result
    for (int i = 0; i<optimalValues.size(); i++) {
        optimalValues[i] /= numValues[i];
    }
    optimal_signal = vnl_vector<double>(&optimalValues[0], (int)optimalValues.size());
    
    return true;
}

/// \brief		Unscented transform of process Sigma points
void UnscentedKalmanFilter::utf(vnl_matrix<double> X, vnl_vector<double> u,
	vnl_vector<double> &y, vnl_matrix<double> &Y, vnl_matrix<double> &P, vnl_matrix<double> &Y1)
{
	// determine number of sigma points
	unsigned int L = X.cols();
	// zero output matrices
	y.fill(0.0); Y.fill(0.0);

	// transform the sigma points and put them as columns in a matrix Y
	for( int k = 0; k < L; k++ )
	{
		vnl_vector<double> xk = X.get_column(k);
		vnl_vector<double> yk(N);
		f(xk,u,yk);
		Y.set_column(k,yk);
		// add each transformed point to the weighted mean
		y = y + Wm.get(0,k)*yk;
	}

	// create a matrix with each column being the weighted mean
	vnl_matrix<double> Ymean(N,L);
	for( int k = 0; k < L; k++ )
		Ymean.set_column(k,y);
	// set the matrix of difference vectors
	Y1 = Y-Ymean;
	// calculate the covariance matrix output
	vnl_matrix<double> WC(L,L,0.0);
	WC.set_diagonal(Wc.get_row(0));
	P = Y1*WC*Y1.transpose();
}

bool rgrsn_ldp::transition_matrix(const vcl_vector<int> & fns,
                                  const vcl_vector<double> & values,
                                  vnl_matrix<double> & transition,
                                  const double resolution)
{
    assert(fns.size() == values.size());
    
    double min_v = *vcl_min_element(values.begin(), values.end());
    double max_v = *vcl_max_element(values.begin(), values.end());
    
    unsigned num_bin = (max_v - min_v)/resolution;
    transition = vnl_matrix<double>(num_bin, num_bin, 0.0);
    vnl_vector<double> column(num_bin, 0.0);
    for (int i = 0; i<fns.size()-1; i++) {
        if (fns[i] + 1 == fns[i+1]) {
            double cur_v = values[i];
            double next_v = values[i+1];
            unsigned cur_bin = value_to_bin_number(min_v, resolution, cur_v, num_bin);
            unsigned next_bin = value_to_bin_number(min_v, resolution, next_v, num_bin);
            transition[next_bin][cur_bin] += 1.0;
            column[cur_bin] += 1.0;
        }
    }
    
    // normalize each column
    for (int r = 0; r < transition.rows(); r++) {
        for (int c = 0; c < transition.cols(); c++) {
            transition[r][c] /= column[c];
        }
    }
    return true;
}

vnl_vector<double> MCLR_SM::Column_Order_Matrix(vnl_matrix<double> mat)
{
		vnl_vector<double> mat_column_ordered;
		mat_column_ordered.set_size(mat.rows()*mat.cols());
		int count = 0;

		for(int j=0;j<mat.cols();++j)	
		{
		for(int i=0;i<mat.rows();++i)
		{
			mat_column_ordered(count) = mat(i,j);
			count++;
		  }
		}
	return mat_column_ordered;
}

/** Applies an affine transform specified by the homogeneous transformation matrix to a vtkPolyData */
void psciob::AffineTransformVTKPolyData( vtkPolyData* inputPolyData, vnl_matrix<double> transformMatrix, vtkPolyData* outputPolyData) {
	if ( transformMatrix.rows() != transformMatrix.cols() ) throw DeformableModelException("transform matrix must be square");
	unsigned int d = transformMatrix.rows()-1;
	vtkSmartPointer<vtkMatrix4x4> matrix = vtkSmartPointer<vtkMatrix4x4>::New();
	vtkSmartPointer<vtkTransform> transform = vtkSmartPointer<vtkTransform>::New();
	vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();

	switch(d) {
	case 2:
		matrix->SetElement(0,0, transformMatrix(0,0)); matrix->SetElement(0,1, transformMatrix(0,1)); matrix->SetElement(0,2, 0); matrix->SetElement(0,3, transformMatrix(0,2)); 
		matrix->SetElement(1,0, transformMatrix(1,0)); matrix->SetElement(1,1, transformMatrix(1,1)); matrix->SetElement(1,2, 0); matrix->SetElement(1,3, transformMatrix(1,2)); 
		matrix->SetElement(2,0,         0           ); matrix->SetElement(2,1,         0           ); matrix->SetElement(2,2, 1); matrix->SetElement(2,3,         0           ); 
		matrix->SetElement(3,0,         0           ); matrix->SetElement(3,1,         0           ); matrix->SetElement(3,2, 0); matrix->SetElement(3,3,         1           ); 
		break;
	case 3:
		for (unsigned i=0 ; i<4 ; i++) {
			for (unsigned j=0 ; j<4 ; j++) { 
				matrix->SetElement(i,j, transformMatrix(i,j)); 
			} 
		}
		break;
	default: throw DeformableModelException("only 2D or 3D transforms are allowed -> 3*3 or 4*4 matrices");
	}

	transform->SetMatrix( matrix );
	transform->Update();

	transformFilter->SetInput( inputPolyData );
	transformFilter->SetTransform( transform );
	transformFilter->Update();

	outputPolyData->ShallowCopy( transformFilter->GetOutput() );
}

void denormalize(vnl_matrix<double>& x, const vnl_vector<double>& centroid, const double scale) {
  int n = x.rows();
  if (n==0) return;
  int d = x.cols();
  for (int i = 0; i < n; ++i) {
    x.set_row(i, x.get_row(i) * scale + centroid);
  }
}

vnl_matrix<double> MCLR_SM::Add_Bias(vnl_matrix<double> data)
{
	vnl_matrix<double> data_with_bias; // Data to be modified

	vnl_vector<double> temp_vector;
	temp_vector.set_size(data.cols());
	temp_vector.fill(1);
	
	data_with_bias.set_size(no_of_features+1,data.cols());
	data_with_bias.set_row(0,temp_vector);
	for(int i =0; i<no_of_features ; ++i)
	{
		data_with_bias.set_row(i+1,data.get_row(i));	
	}
		
	return data_with_bias;
}

double GaussTransform(const vnl_matrix<double>& A,
    const vnl_matrix<double>& B, double scale,
    vnl_matrix<double>& gradient) {
  // assert A.cols() == B.cols()
  return GaussTransform(A.data_block(), B.data_block(),
      A.rows(), B.rows(), A.cols(), scale,
      gradient.data_block());
}

void normalize_same(vnl_matrix<double>& x,
    vnl_vector<double>& centroid, double& scale) {
  int n = x.rows();
  if (n==0) return;
  int d = x.cols();
  for (int i = 0; i < n; ++i) {
    x.set_row(i, (x.get_row(i) - centroid) / scale);
  }
}

int select_points(const vnl_matrix<double>&pts,
    const std::vector<int>&index, vnl_matrix<double>& selected) {
  int n = index.size();
  int d = pts.cols();
  selected.set_size(n,d);
  for (int i = 0; i < n; ++i) {
    selected.update(pts.extract(1, d, index[i]), i);
  }
  return n;
}