C++ Matrix2D类代码示例

void AllSubsets::reset(Matrix2D* m, Matrix2D& points, Vector& mu)
{
  reset();

  //create first subset matrix
  for(int j = 0; j <= sizeOfSubsample; j++)
    m->setValue(0, j, 1.0);

  for(int i = 1; i <= sizeOfSubsample; i++)
  {
    m->setValue(i, 0, mu.getValue(i-1));
  }

  for(int j = 0; j < sizeOfSubsample - 1; j++)
  {
    selectedElements[j] = j;
    for(int i = 0; i < sizeOfSubsample; i++)
    {
      m->setValue(i+1, j+1, points.getValue(i,j));
    }
  }

  selectedElements[numberOfElements - 1] = sizeOfSubsample - 1;
  for(int i = 0; i < sizeOfSubsample; i++)
    m->setValue(i+1, sizeOfSubsample, points.getValue(i, numberOfElements - 1));
}

//#define DEBUG
void ProgSSNR::run()
{
    show();
    produceSideInfo();

    Matrix2D<double> output;
    if (!radial_avg)
    {
        if (!generate_VSSNR)
            estimateSSNR(1, output);
        else
            estimateSSNR(2, output);
        if (fn_out == "")
        	fn_out=fn_S.insertBeforeExtension("_SSNR").removeLastExtension().addExtension("xmd");
    }
    else
    {
        radialAverage(output);
        if (fn_out == "")
        	fn_out=fn_VSSNR.insertBeforeExtension("_radial_avg").removeLastExtension().addExtension("xmd");
    }
#ifdef DEBUG
    output.write(fn_out);
#endif
    MetaData MD;
    for (size_t i=1; i<MAT_YSIZE(output); ++i)
    {
    	size_t id=MD.addObject();
    	MD.setValue(MDL_RESOLUTION_FREQ,output(i,1),id);
  		MD.setValue(MDL_RESOLUTION_SSNR,output(i,2),id);
    	MD.setValue(MDL_RESOLUTION_FREQREAL,1.0/output(i,1),id);
    }
    MD.write(fn_out);
}

int MatrixBinder::setTy(lua_State* L)
{
	Binder binder(L);
	Matrix2D* matrix = static_cast<Matrix2D*>(binder.getInstance("Matrix", 1));

	matrix->setTy(luaL_checknumber(L, 2));

	return 0;
}

void TextureTest::ProcessEvent(SDL_Event* event)
{
    switch(event->type)
    {
        case SDL_MOUSEBUTTONDOWN:
        {
            if(event->button.button == 4)
            {

            }
            else if(event->button.button == 5)
            {

            }
            break;
        }
        case SDL_KEYDOWN:
        {
            if(event->key.keysym.sym == SDLK_UP)
            {
                Vector2D oldSteer = v->GetVelocity();
                Vector2D temp = v->GetVelocity();
                temp.Normalize();
                v->SetVelocity(oldSteer+temp*16);
                //v->SetHeading(v->GetVelocity());
            }
            if(event->key.keysym.sym == SDLK_DOWN)
            {
                Vector2D oldSteer = v->GetVelocity();
                Vector2D temp = v->GetVelocity();
                temp.Normalize();
                v->SetVelocity(oldSteer-temp*163);
            }
            if(event->key.keysym.sym == SDLK_LEFT)
            {
                Matrix2D mat;
                //Vector2D vec = v->GetVelocity();
                Vector2D vecH = v->GetVelocity();
                //Vec2DRotateAroundO(vecH,angle);
                mat.Rotate(angle);
//
               //mat.Rotate(angle);
                mat.TransformVector(vecH);
//              / mat.TransformVector(vecH);
//                //v->SetVelocity(vec);
//                v->RotateHeadingToFacePosition(vecH);
                //v->SetSteeringForce(vecH);
                v->SetVelocity(vecH);
                angle=0.07;

            }
            break;
        }
        //case SDL_KEYUP
    }
}

void
GaussJordan::print(Matrix2D const & A) const {
    std::cout << std::endl;
    for (IMatrix2D::size_type i = 0; i < A.rows(); ++i) {
        IMatrix2D::size_type row = logicalToPhysicalRowIndex(i);
        for (IMatrix2D::size_type j = 0; j < A.cols(); ++j) {
            std::cout << std::setw(10) << A(row, j);
        }
        std::cout << std::endl;
    }
}

/* Interface to numerical recipes: svbksb ---------------------------------- */
void svbksb(Matrix2D<double> &u, Matrix1D<double> &w, Matrix2D<double> &v,
            Matrix1D<double> &b, Matrix1D<double> &x)
{
    // Call to the numerical recipes routine. Results will be stored in X
    svbksb(u.adaptForNumericalRecipes2(),
           w.adaptForNumericalRecipes(),
           v.adaptForNumericalRecipes2(),
           u.mdimy, u.mdimx,
           b.adaptForNumericalRecipes(),
           x.adaptForNumericalRecipes());
}

Matrix2D GetImageLineMatrix(double x1, double y1, double x2, double y2, const Image& image)
{
	Matrix2D m;
	Size sz = image.GetSize();
	m.scale(agg::calc_distance(x1, y1, x2, y2) / sz.cx);
	if(fabs(x2 - x1) < fabs(y2 - y1) * 1e-6)
		m.rotate(y2 > y1 ? M_PI_2 : -M_PI_2);
	else
		m.rotate(atan((y2 - y1) / (x2 - x1)));
	m.translate(x1, y1);
	return m;
}

// Get matrix ==============================================================
void SymList::get_matrices(int i, Matrix2D<DOUBLE> &L, Matrix2D<DOUBLE> &R)
const
{
    int k, l;
    L.initZeros(4, 4);
    R.initZeros(4, 4);
    for (k = 4 * i; k < 4*i + 4; k++)
        for (l = 0; l < 4; l++)
        {
            L(k - 4*i, l) = __L(k, l);
            R(k - 4*i, l) = __R(k, l);
        }
}

	void optimiseTransformationMatrixContinuous()
	{
		// Get coordinates of all pairs:
		Matrix2D<double> Au, Bt;
	    Au.initZeros(3, 3);
	    Bt.initZeros(3, 3);
	    Pass.initZeros(4,4);

	    // Add all pairs to dependent matrices (adapted from add_point in Xmipps micrograph_mark main_widget_mark.cpp)
		for (int t = 0; t < pairs_t2u.size(); t++)
		{
			int u = pairs_t2u[t];
			if (u >= 0)
	        {
				Au(0, 0) += (double)(p_unt[2*u] * p_unt[2*u]);
				Au(0, 1) += (double)(p_unt[2*u] * p_unt[2*u+1]);
				Au(0, 2) += (double)(p_unt[2*u]);
				Au(1, 0) = Au(0, 1);
				Au(1, 1) += (double)(p_unt[2*u+1] * p_unt[2*u+1]);
				Au(1, 2) += (double)(p_unt[2*u+1]);
				Au(2, 0) = Au(0, 2);
				Au(2, 1) = Au(1, 2);
				Au(2, 2) += 1.;

				Bt(0, 0) += (double)(p_til[2*t] * p_unt[2*u]);
				Bt(0, 1) += (double)(p_til[2*t+1] * p_unt[2*u]);
				Bt(0, 2) = Au(0, 2);
				Bt(1, 0) += (double)(p_til[2*t] * p_unt[2*u+1]);
				Bt(1, 1) += (double)(p_til[2*t+1] * p_unt[2*u+1]);
				Bt(1, 2) = Au(1, 2);
				Bt(2, 0) += (double)(p_til[2*t]);
				Bt(2, 1) += (double)(p_til[2*t+1]);
				Bt(2,2) += 1.;
	        }
	    }

	    // Solve equations
	    solve(Au, Bt, Pass);
	    Pass = Pass.transpose();
	    std::cout << " Optimised passing matrix= " << Pass << std::endl;
	    //These values can be complete CRAP. Better not show them at all....
	    //double rotp, tiltp, psip;
	    //tiltp = acos(Pass(1,1));
	    //rotp = acos(Pass(1,0)/sin(tiltp));
	    //psip = acos(Pass(0,1)/-sin(tiltp));
	    //std::cout << " Optimised tilt angle= " << RAD2DEG(tiltp) << std::endl;
	    //std::cout << " Optimised in-plane rot angles= " << RAD2DEG(rotp) <<" and "<< RAD2DEG(psip) << std::endl;
	    // Map using the new matrix
	    mapOntoTilt();

	}

double JPetRecoImageTools::calculateProjection(const Matrix2D& emissionMatrix, double angle, int scanNumber, int nScans,
    InterpolationFunc& interpolationFunction)
{
  int N = scanNumber - nScans / 2 ;
  const int kInputMatrixSize = emissionMatrix.size();
  //if no. nScans is greater than the image width, then scale will be <1
  const double scale = kInputMatrixSize / nScans;
  const double kSin45or125deg = std::sqrt(2) / 2; /// sin(45) deg
  const double kEpsilon = 0.0000001;
  const double kDegToRad = M_PI / 180.;

  double sin = std::sin(angle * kDegToRad - M_PI / 2.);
  sin = setToZeroIfSmall(sin, kEpsilon);
  double cos = std::cos(angle * kDegToRad - M_PI / 2.);
  cos = setToZeroIfSmall(cos, kEpsilon);

  double a = 0.;
  double b = 0.;
  /// The line over which we integrate is perpendicular to any line with the slope = tg(angle), so it is always  -1/tg(angle).
  /// If the angle is between (45 to 125)
  /// we use  y = a * x + b, and we iterate over rows of the matrix (x).
  /// If the angle is between [0 to 45 ] or [125 to 180]
  /// we use x = a* y +b, and we iterate over columns of the matrix (y)
  bool angleRange45To125 = std::abs(sin) > kSin45or125deg;
  double divided = 1.;
  std::function<double(int, int)> matrixGet;

  if (angleRange45To125) {
    assert(sin);
    a = -cos / sin;
    b = (N - cos - sin) / sin;
    b *= scale;
    matrixGet  = matrixGetterFactory(emissionMatrix, false); // The matrix  elements will be taken as (x,y).
    divided = std::abs(sin);
  } else {
    assert(cos);
    a = -sin / cos;
    b = (N - cos - sin) / cos;
    b *= scale;
    matrixGet  = matrixGetterFactory(emissionMatrix, true); // The matrix  elements will be taken as (y, x) - transposed.
    divided = std::abs(cos);
  }
  const int kMatrixCenter = emissionMatrix.size() / 2;
  double value = 0.;
  for (auto i = -kMatrixCenter; i < kMatrixCenter; i++) {
    value += interpolationFunction(i + kMatrixCenter , a * i + b + kMatrixCenter, matrixGet);
  }
  value /= divided;
  return value;
}

void
GaussJordan::rearrangeDueToPivoting(Matrix2D & A, Matrix2D & AInverse, Vector & rhs) const {
    decltype(partial_pivoting_map_) physical_map(partial_pivoting_map_.size());
    for (auto i = 0; i < partial_pivoting_map_.size(); ++i) {
        physical_map[i] = physicalToLogicalRowIndex(i);
    }

    for (auto i = 0; i < physical_map.size(); ++i) {
        auto j = std::distance(std::begin(physical_map), std::find(std::begin(physical_map), std::end(physical_map), i));
        if (i == j)
            continue;

        // Swap rows
        for (auto col = 0; col < A.cols(); ++col) {
            std::swap(AInverse(i, col), AInverse(j, col));
            std::swap(A(i, col), A(j, col));
        }

//         AInverse.print();
//         print(AInverse);

        // Swap rows on the r.h.s.
        std::swap(rhs(i), rhs(j));
        std::swap(physical_map[i], physical_map[j]);
    }

//     AInverse.print();
//     print(AInverse);
}

JPetRecoImageTools::Matrix2DProj JPetRecoImageTools::sinogram(Matrix2D& emissionMatrix,
    int nViews, int nScans,
    double angleBeg, double angleEnd,
    InterpolationFunc interpolationFunction,
    RescaleFunc rescaleFunc,
    int minCutoff,
    int scaleFactor
                                                             )
{
  assert(emissionMatrix.size() > 0);
  assert(emissionMatrix.size() == emissionMatrix[0].size());
  assert(nViews > 0);
  assert(nScans > 0);
  assert(angleBeg < angleEnd);
  assert(minCutoff < scaleFactor);

  //create vector of size nViews, initialize it with vector of size nScans
  Matrix2DProj proj(nViews, std::vector<double>(nScans));

  float stepsize = (angleEnd - angleBeg) / nViews;
  assert(stepsize > 0); //maybe != 0 ?

  int viewIndex = 0;
  for (auto phi = angleBeg; phi < angleEnd; phi = phi + stepsize, viewIndex++) {
    for (auto scanNumber = 0; scanNumber < nScans; scanNumber++) {
      proj[viewIndex][nScans  - 1 - scanNumber] = JPetRecoImageTools::calculateProjection(emissionMatrix,
          phi, scanNumber, nScans,
          interpolationFunction);
    }
  }
  rescaleFunc(proj, minCutoff, scaleFactor);
  return proj;
}

Matrix2D operator *( const Matrix2D &first, const Matrix2D &second ) {
   Matrix2D returnVal;
   
   for( int n = 0; n < Matrix2D::SIDE_LENGTH; n++ ) {
      for( int p = 0; p < Matrix2D::SIDE_LENGTH; p++ ) {
         float val = 0.0;
         
         for( int m = 0; m < Matrix2D::SIDE_LENGTH; m++ ) {
            val += first.Get( n, m ) + second.Get( m, p );
         }
         
         returnVal.Set( n, p, val );
      }
   }
   
   return returnVal;
}

std::function<double(int, int)> JPetRecoImageTools::matrixGetterFactory(const Matrix2D& emissionMatrix, bool isTransposed)
{
  if (!isTransposed) {
    return [& emissionMatrix](int i, int j) {
      if (i >= 0 && i < (int) emissionMatrix[0].size() && j >= 0 && j < (int) emissionMatrix.size() ) {
        return emissionMatrix[i][j];
      } else {
        return 0;
      }
    };
  } else {
    return [& emissionMatrix](int i, int j) {
      if (i >= 0 && i < (int) emissionMatrix.size() && j >= 0 && j < (int) emissionMatrix[0].size() ) {
        return emissionMatrix[j][i];
      } else {
        return 0;
      }
    };
  }
}

double FieldVariableDescriptor::GetComponent(const Matrix2D & m) const
{
	//EK 2012-03-02 assert bug fix - if component_index_1 == FVCI_magnitude, component_index_2 does not matter
	assert(component_index_1 != FVCI_none && component_index_2 != FVCI_magnitude);
	assert(component_index_2 != FVCI_none || component_index_1 == FVCI_magnitude);
	if(component_index_1 == FVCI_magnitude)
		return sqrt(m.InnerProduct(m));
	if(component_index_1 == FVCI_z || component_index_2 == FVCI_z)
		return FIELD_VARIABLE_NO_VALUE;		// this may happen if there are 3D and 2D elements mixed in the system
	return m(component_index_1, component_index_2);
}