import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes a basis (the principal components) from the most dominant eigenvectors.
 *
 * @param numComponents Number of vectors it will use to describe the data.  Typically much
 *                      smaller than the number of elements in the input vector.
 */
public void computeBasis(int numComponents) {
    Logger.getAnonymousLogger().log(Level.INFO, "Compute basis to get principal components.");

    validateData(numComponents);
    setNumPrincipalComponents(numComponents);
    computeNormalizedMean();

    // Compute SVD and save time by not computing U
    SingularValueDecomposition<DenseMatrix64F> svd =
            DecompositionFactory.svd(getDataMatrix().numRows, getDataMatrix().numCols, false, true, false);
    if (!svd.decompose(getDataMatrix())) {
        throw new IllegalStateException("SVD failure. Can't compute principal components!");
    }

    setPrincipalComponentSubspace(svd.getV(null, true));
    final DenseMatrix64F singularDiagonalMatrix = svd.getW(null);

    // Singular values are in an arbitrary order initially. We ask for principal components subspace to be transposed.
    SingularOps.descendingOrder(null, false, singularDiagonalMatrix, getPrincipalComponentSubspace(), true);

    // strip off unneeded components and find the basis
    getPrincipalComponentSubspace().reshape(numComponents, getMean().length, true);
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes the SVD of A and extracts the homography matrix from its null space
 */
protected boolean computeH(DenseMatrix64F A, DenseMatrix64F H) {
	if( !svd.decompose(A) )
		return true;

	if( A.numRows > 8 )
		SingularOps.nullVector(svd,true,H);
	else {
		// handle a special case since the matrix only has 8 singular values and won't select
		// the correct column
		DenseMatrix64F V = svd.getV(null,false);
		SpecializedOps.subvector(V, 0, 8, V.numCols, false, 0, H);
	}

	return false;
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes the null space of A and extracts the transform.
 */
private void computeTransform( DenseMatrix64F A ) {
	if( !svd.decompose(A) )
		throw new RuntimeException("SVD failed?");

	SingularOps.nullVector(svd,true,x);

	DenseMatrix64F R = motion.getR();
	Vector3D_F64 T = motion.getT();

	// extract the results
	System.arraycopy(x.data,0,R.data,0,9);

	T.x = x.data[9];
	T.y = x.data[10];
	T.z = x.data[11];
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes the decomposition from an essential matrix.
 *
 * @param E essential matrix
 */
public void decompose( DenseMatrix64F E ) {
	if( svd.inputModified() ) {
		E_copy.set(E);
		E = E_copy;
	}

	if( !svd.decompose(E))
		throw new RuntimeException("Svd some how failed");

	U = svd.getU(U,false);
	V = svd.getV(V,false);
	S = svd.getW(S);

	SingularOps.descendingOrder(U,false,S,V,false);

	decompose(U, S, V);
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Projects the found estimate of E onto essential space.
 *
 * @return true if svd returned true.
 */
protected boolean projectOntoEssential( DenseMatrix64F E ) {
	if( !svdConstraints.decompose(E) ) {
		return false;
	}
	svdV = svdConstraints.getV(svdV,false);
	svdU = svdConstraints.getU(svdU,false);
	svdS = svdConstraints.getW(svdS);

	SingularOps.descendingOrder(svdU, false, svdS, svdV, false);

	// project it into essential space
	// the scale factor is arbitrary, but the first two singular values need
	// to be the same.  so just set them to one
	svdS.unsafe_set(0, 0, 1);
	svdS.unsafe_set(1, 1, 1);
	svdS.unsafe_set(2, 2, 0);

	// recompute F
	CommonOps.mult(svdU, svdS, temp0);
	CommonOps.multTransB(temp0,svdV, E);

	return true;
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Projects the found estimate of F onto Fundamental space.
 *
 * @return true if svd returned true.
 */
protected boolean projectOntoFundamentalSpace( DenseMatrix64F F ) {
	if( !svdConstraints.decompose(F) ) {
		return false;
	}
	svdV = svdConstraints.getV(svdV,false);
	svdU = svdConstraints.getU(svdU,false);
	svdS = svdConstraints.getW(svdS);

	SingularOps.descendingOrder(svdU, false, svdS, svdV, false);

	// the smallest singular value needs to be set to zero, unlike
	svdS.set(2, 2, 0);

	// recompute F
	CommonOps.mult(svdU, svdS, temp0);
	CommonOps.multTransB(temp0,svdV, F);

	return true;
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes the SVD of A and extracts the essential/fundamental matrix from its null space
 */
protected boolean process(DenseMatrix64F A, DenseMatrix64F F ) {
	if( !svdNull.decompose(A) )
		return true;

	if( A.numRows > 8 )
		SingularOps.nullVector(svdNull,true,F);
	else {
		// handle a special case since the matrix only has 8 singular values and won't select
		// the correct column
		DenseMatrix64F V = svdNull.getV(null,false);
		SpecializedOps.subvector(V, 0, 8, V.numCols, false, 0, F);
	}

	return false;
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Given a set of homographies computed from a sequence of images that observe the same
 * plane it estimates the camera's calibration.
 *
 * @param homographies Homographies computed from observations of the calibration grid.
 */
public void process( List<DenseMatrix64F> homographies ) {
	if( assumeZeroSkew ) {
		if( homographies.size() < 2 )
			throw new IllegalArgumentException("At least two homographies are required");
	} else if( homographies.size() < 3 ) {
		throw new IllegalArgumentException("At least three homographies are required");
	}

	if( assumeZeroSkew ) {
		setupA_NoSkew(homographies);
		if( !svd.decompose(A) )
			throw new RuntimeException("SVD failed");
		if( homographies.size() == 2 ) {
			DenseMatrix64F V = svd.getV(null,false);
			SpecializedOps.subvector(V, 0, 4, V.numRows, false, 0, b);
		} else {
			SingularOps.nullVector(svd,true,b);
		}
		computeParam_ZeroSkew();
	} else {
		setupA(homographies);
		if( !svd.decompose(A) )
			throw new RuntimeException("SVD failed");
		SingularOps.nullVector(svd,true,b);
		computeParam();
	}
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * <p>
 * Given N observations of the same point from two views and a known motion between the
 * two views, triangulate the point's position in camera 'b' reference frame.
 * </p>
 * <p>
 * Modification of [1] to be less generic and use calibrated cameras.
 * </p>
 *
 * @param observations Observation in each view in normalized coordinates. Not modified.
 * @param worldToView Transformations from world to the view.  Not modified.
 * @param found Output, the found 3D position of the point.  Modified.
 */
public void triangulate( List<Point2D_F64> observations ,
						 List<Se3_F64> worldToView ,
						 Point3D_F64 found ) {
	if( observations.size() != worldToView.size() )
		throw new IllegalArgumentException("Number of observations must match the number of motions");
	
	final int N = worldToView.size();
	
	A.reshape(2*N,4,false);
	
	int index = 0;
	
	for( int i = 0; i < N; i++ ) {
		index = addView(worldToView.get(i),observations.get(i),index);
	}

	if( !svd.decompose(A) )
		throw new RuntimeException("SVD failed!?!?");

	SingularOps.nullVector(svd,true,v);

	double w = v.get(3);
	found.x = v.get(0)/w;
	found.y = v.get(1)/w;
	found.z = v.get(2)/w;
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
	 * Computes M'*M and finds the null space.  The 4 eigenvectors with the smallest eigenvalues are found
	 * and the null points extracted from them.
	 */
	protected void extractNullPoints( DenseMatrix64F M )
	{
		// compute MM and find its null space
		MM.reshape(M.numRows,M.numRows,false);
		CommonOps.multTransB(M, M, MM);

		if( !svd.decompose(MM) )
			throw new IllegalArgumentException("SVD failed?!?!");

		double []singularValues = svd.getSingularValues();
		DenseMatrix64F V = svd.getV(null,false);

		SingularOps.descendingOrder(null,false,singularValues,3,V,false);

		// extract null points from the null space
		for( int i = 0; i < numControl; i++ ) {
			int column = M.numRows-1-i;
//			nullPts[i].clear();
			for( int j = 0; j < numControl; j++ ) {
				Point3D_F64 p = nullPts[i].get(j);
//				Point3D_F64 p = new Point3D_F64();
				p.x = V.get(j*3+0,column);
				p.y = V.get(j*3+1,column);
				p.z = V.get(j*3+2,column);
//				nullPts[i].add(p);
			}
		}
	}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes a trifocal tensor which minimizes the algebraic error given the
 * two epipoles and the linear constraint matrix.  The epipoles are from a previously
 * computed trifocal tensor.
 *
 * @param e2 Epipole of first image in the second image
 * @param e3 Epipole of first image in the third image
 * @param A Linear constraint matrix for trifocal tensor created from image observations.
 */
public void process( Point3D_F64 e2 , Point3D_F64 e3 , DenseMatrix64F A ) {
	// construct the linear system that the solution which solves the unknown square
	// matrices in the camera matrices
	constructE(e2, e3);

	// Computes U, which is used to map the 18 unknowns onto the 27 trifocal unknowns
	svdU.decompose(E);
	svdU.getU(U, false);

	// Copy the parts of U which correspond to the non singular parts if the SVD
	// since there are only really 18-nullity unknowns due to linear dependencies
	SingularOps.descendingOrder(U,false,svdU.getSingularValues(),svdU.numberOfSingularValues(),null,false);
	int rank = SingularOps.rank(svdU, 1e-13);
	Up.reshape(U.numRows,rank);
	CommonOps.extract(U,0,U.numRows,0,Up.numCols,Up,0,0);

	// project the linear constraint matrix into this subspace
	AU.reshape(A.numRows,Up.numCols);
	CommonOps.mult(A,Up,AU);

	// Extract the solution of ||A*U*x|| = 0 from the null space
	svdV.decompose(AU);

	xp.reshape(rank,1);
	SingularOps.nullVector(svdV,true,xp);

	// Translate the solution from the subspace and into a valid trifocal tensor
	CommonOps.mult(Up,xp,vectorT);

	// the sign of vectorT is arbitrary, but make it positive for consistency
	if( vectorT.data[0] > 0 )
		CommonOps.changeSign(vectorT);
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes the null space of the linear system to find the trifocal tensor
 */
protected boolean solveLinearSystem() {
	if( !svdNull.decompose(A) )
		return false;

	SingularOps.nullVector(svdNull,true,vectorizedSolution);

	solutionN.convertFrom(vectorizedSolution);

	return true;
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Extracts the epipoles from the trifocal tensor.  Extracted epipoles will have a norm of 1
 * as an artifact of using SVD.
 *
 * @param tensor Input: Trifocal tensor.  Not Modified
 * @param e2  Output: Epipole in image 2. Homogeneous coordinates. Modified
 * @param e3  Output: Epipole in image 3. Homogeneous coordinates. Modified
 */
public void process( TrifocalTensor tensor , Point3D_F64 e2 , Point3D_F64 e3 ) {
	svd.decompose(tensor.T1);
	SingularOps.nullVector(svd, true, v1);
	SingularOps.nullVector(svd, false,u1);

	svd.decompose(tensor.T2);
	SingularOps.nullVector(svd,true,v2);
	SingularOps.nullVector(svd,false,u2);

	svd.decompose(tensor.T3);
	SingularOps.nullVector(svd,true,v3);
	SingularOps.nullVector(svd,false,u3);

	for( int i = 0; i < 3; i++ ) {
		U.set(i,0,u1.get(i));
		U.set(i,1,u2.get(i));
		U.set(i,2,u3.get(i));

		V.set(i, 0, v1.get(i));
		V.set(i, 1, v2.get(i));
		V.set(i, 2, v3.get(i));
	}

	svd.decompose(U);
	SingularOps.nullVector(svd, false, tempE);
	e2.set(tempE.get(0), tempE.get(1), tempE.get(2));

	svd.decompose(V);
	SingularOps.nullVector(svd, false, tempE);
	e3.set(tempE.get(0), tempE.get(1), tempE.get(2));
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes a basis (the principal components) from the most dominant eigenvectors.
 *
 * @param numComponents Number of vectors it will use to describe the data.  Typically much
 * smaller than the number of elements in the input vector.
 */
public void computeBasis( int numComponents ) {
    if( numComponents > A.getNumCols() )
        throw new IllegalArgumentException("More components requested that the data's length.");
    if( sampleIndex != A.getNumRows() )
        throw new IllegalArgumentException("Not all the data has been added");
    if( numComponents > sampleIndex )
        throw new IllegalArgumentException("More data needed to compute the desired number of components");

    this.numComponents = numComponents;

    // compute the mean of all the samples
    for( int i = 0; i < A.getNumRows(); i++ ) {
        for( int j = 0; j < mean.length; j++ ) {
            mean[j] += A.get(i,j);
        }
    }
    for( int j = 0; j < mean.length; j++ ) {
        mean[j] /= A.getNumRows();
    }

    // subtract the mean from the original data
    for( int i = 0; i < A.getNumRows(); i++ ) {
        for( int j = 0; j < mean.length; j++ ) {
            A.set(i,j,A.get(i,j)-mean[j]);
        }
    }

    // Compute SVD and save time by not computing U
    SingularValueDecomposition<DenseMatrix64F> svd =
            DecompositionFactory.svd(A.numRows, A.numCols, false, true, false);
    if( !svd.decompose(A) )
        throw new RuntimeException("SVD failed");

    V_t = svd.getV(null,true);
    DenseMatrix64F W = svd.getW(null);

    // Singular values are in an arbitrary order initially
    SingularOps.descendingOrder(null,false,W,V_t,true);

    // strip off unneeded components and find the basis
    V_t.reshape(numComponents,mean.length,true);
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes a basis (the principle components) from the most dominant eigenvectors.
 *
 * @param numComponents Number of vectors it will use to describe the data.  Typically much
 * smaller than the number of elements in the input vector.
 */
public void computeBasis( int numComponents ) {
   if( numComponents > A.getNumCols() )
      throw new IllegalArgumentException("More components requested that the data's length.");
   if( sampleIndex != A.getNumRows() )
      throw new IllegalArgumentException("Not all the data has been added");
   if( numComponents > sampleIndex )
      throw new IllegalArgumentException("More data needed to compute the desired number of components");

   this.numComponents = numComponents;

   // compute the mean of all the samples
   for( int i = 0; i < A.getNumRows(); i++ ) {
      for( int j = 0; j < mean.length; j++ ) {
         mean[j] += A.get(i,j);
      }
   }
   for( int j = 0; j < mean.length; j++ ) {
      mean[j] /= A.getNumRows();
   }

   // subtract the mean from the original data
   for( int i = 0; i < A.getNumRows(); i++ ) {
      for( int j = 0; j < mean.length; j++ ) {
         A.set(i,j,A.get(i,j)-mean[j]);
      }
   }

   // Compute SVD and save time by not computing U
   SingularValueDecomposition<DenseMatrix64F> svd =
          DecompositionFactory.svd(A.numRows, A.numCols, false, true, false);
   if( !svd.decompose(A) )
      throw new RuntimeException("SVD failed");

   V_t = svd.getV(null,true);
   DenseMatrix64F W = svd.getW(null);

   // Singular values are in an arbitrary order initially
   SingularOps.descendingOrder(null,false,W,V_t,true);

   // strip off unneeded components and find the basis
   V_t.reshape(numComponents,mean.length,true);
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Selects control points along the data's axis and the data's centroid.  If the data is determined
 * to be planar then only 3 control points are selected.
 *
 * The data's axis is determined by computing the covariance matrix then performing SVD.  The axis
 * is contained along the
 */
public void selectWorldControlPoints(List<Point3D_F64> worldPts, FastQueue<Point3D_F64> controlWorldPts) {

	meanWorldPts = UtilPoint3D_F64.mean(worldPts);

	// covariance matrix elements, summed up here for speed
	double c11=0,c12=0,c13=0,c22=0,c23=0,c33=0;

	final int N = worldPts.size();
	for( int i = 0; i < N; i++ ) {
		Point3D_F64 p = worldPts.get(i);

		double dx = p.x- meanWorldPts.x;
		double dy = p.y- meanWorldPts.y;
		double dz = p.z- meanWorldPts.z;

		c11 += dx*dx;c12 += dx*dy;c13 += dx*dz;
		c22 += dy*dy;c23 += dy*dz;
		c33 += dz*dz;
	}
	c11/=N;c12/=N;c13/=N;c22/=N;c23/=N;c33/=N;

	DenseMatrix64F covar = new DenseMatrix64F(3,3,true,c11,c12,c13,c12,c22,c23,c13,c23,c33);

	// find the data's orientation and check to see if it is planar
	svd.decompose(covar);
	double []singularValues = svd.getSingularValues();
	DenseMatrix64F V = svd.getV(null,false);

	SingularOps.descendingOrder(null,false,singularValues,3,V,false);

	// planar check
	if( singularValues[0]<singularValues[2]*1e13 ) {
		numControl = 4;
	} else {
		numControl = 3;
	}

	// put control points along the data's major axises
	controlWorldPts.reset();
	for( int i = 0; i < numControl-1; i++ ) {
		double m = Math.sqrt(singularValues[1])*magicNumber;

		double vx = V.unsafe_get(0, i)*m;
		double vy = V.unsafe_get(1, i)*m;
		double vz = V.unsafe_get(2, i)*m;

		controlWorldPts.grow().set(meanWorldPts.x + vx, meanWorldPts.y + vy, meanWorldPts.z + vz);
	}
	// set a control point to be the centroid
	controlWorldPts.grow().set(meanWorldPts.x, meanWorldPts.y, meanWorldPts.z);
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Decomposed the provided homography matrix into its R,T/d,N components. Four
 * solutions will be produced and can be accessed with {@link #getSolutionsN()} and
 * {@link #getSolutionsSE()}.
 *
 * @param H Homography matrix.  Not modified.
 */
public void decompose( DenseMatrix64F H ) {
	if( !svd.decompose(H) )
		throw new RuntimeException("SVD failed somehow");

	DenseMatrix64F V = svd.getV(null,false);
	DenseMatrix64F S = svd.getW(null);

	SingularOps.descendingOrder(null,false,S, V,false);

	// COMMENT: The smallest singular value should be zero so I'm not sure why
	// that is not assumed here.  Seen the same strategy done in a few papers
	// Maybe that really isn't always the case?
	double s0 = S.get(0,0)*S.get(0,0);
	// the middle singular value is known to be one
	double s2 = S.get(2,2)*S.get(2,2);

	v2.set(V.get(0,1),V.get(1,1),V.get(2,1));

	double a = Math.sqrt(1-s2);
	double b = Math.sqrt(s0-1);
	double div = Math.sqrt(s0-s2);

	for( int i = 0; i < 3; i++ ) {
		double e1 = (a*V.get(i,0) + b*V.get(i,2))/div;
		double e2 = (a*V.get(i,0) - b*V.get(i,2))/div;

		u1.setIndex(i,e1);
		u2.setIndex(i,e2);
	}

	setU(U1, v2, u1);
	setU(U2, v2, u2);
	setW(W1, H , v2, u1);
	setW(W2, H , v2, u2);

	// create the four solutions
	createSolution(W1, U1, u1,H,solutionsSE.get(0), solutionsN.get(0));
	createSolution(W2, U2, u2,H,solutionsSE.get(1), solutionsN.get(1));
	createMirrorSolution(0,2);
	createMirrorSolution(1,3);
}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 * Computes a basis (the principle components) from the most dominant eigenvectors.
 */
public void computeBasis() {
	if (sampleIndex != numSamples)
		throw new IllegalArgumentException("Not all the data has been added");
	if (numComponents > numSamples)
		throw new IllegalArgumentException(
				"More data needed to compute the desired number of components");

	means = new DenseMatrix64F(sampleSize, 1);
	// compute the mean of all the samples
	for (int i = 0; i < numSamples; i++) {
		for (int j = 0; j < sampleSize; j++) {
			double val = means.get(j);
			means.set(j, val + A.get(i, j));
		}
	}
	for (int j = 0; j < sampleSize; j++) {
		double avg = means.get(j) / numSamples;
		means.set(j, avg);
	}

	// subtract the mean from the original data
	for (int i = 0; i < numSamples; i++) {
		for (int j = 0; j < sampleSize; j++) {
			A.set(i, j, A.get(i, j) - means.get(j));
		}
	}

	// compute SVD and save time by not computing U
	SingularValueDecomposition<DenseMatrix64F> svd = DecompositionFactory.svd(numSamples, sampleSize,
			false, true, compact);
	if (!svd.decompose(A))
		throw new RuntimeException("SVD failed");

	V_t = svd.getV(null, true);
	W = svd.getW(null);

	// singular values are in an arbitrary order initially and need to be sorted in descending order
	SingularOps.descendingOrder(null, false, W, V_t, true);

	// strip off unneeded components and find the basis
	V_t.reshape(numComponents, sampleSize, true);

}
 

import org.ejml.ops.SingularOps; //导入依赖的package包/类
/**
 *
 * @param dataMatrix Set of all data vectors, one vector per column
 * @throws Exception
 */
public PCAviaSVD(double[][] dataMatrix) throws Exception {//dataMatrix=[x1|x2|..|xN], xi  D-dim
    int D=dataMatrix.length;
    int N=dataMatrix[0].length;
    
    //evaluating means
    means = new double[D];
    System.out.print("PCA is evaluating mean...");
    for (int i = 0; i < N; i++)
        for (int j = 0; j < D; j++)
            means[j] += dataMatrix[j][i];
    
    for (int i = 0; i < D; i++)
        means[i] /= N;
    
    System.out.println("done");
    //data centering
    DenseMatrix64F X= new DenseMatrix64F(D,N);
    for(int i1=0; i1<D;i1++)
        for(int i2=0;i2<N; i2++)
            X.set(i1, i2, dataMatrix[i1][i2]-means[i1]);
    
    //evaluating SVD
    System.out.print("PCA is evaluating principal components via SVD...");
    SingularValueDecomposition<DenseMatrix64F> svd=DecompositionFactory.svd(D, N, true, false, false);//compute only U
    
    
    long start = System.currentTimeMillis();
    if(!svd.decompose(X))
        throw new Exception("SVD failed");
    
    System.out.println("Time: " + (System.currentTimeMillis()-start)/1000 );
    
    DenseMatrix64F U=svd.getU(null, true);//U trasposed
    DenseMatrix64F W= svd.getW(null);
    
    //sorting principal components
    SingularOps.descendingOrder(U,true , W, null, false);
    
    principalComponentsMatrix=U;
    System.out.println("done");
}