import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
/**
 * Get the closest annotation to a given display point.
 *
 * @param x
 * @param y
 * @param maxDistance Max distance to consider for the search. If no annotation exists in the search area, null is returned.
 * @return
 */
private Annotation getClosestAnnotation(final int x, final int y, final double maxDistance) {

	final RealPoint pos = new RealPoint(3);
	viewer.displayToGlobalCoordinates(x, y, pos);
	System.out.println("clicked at global coordinates " + pos);

	final List< Annotation > closest = annotations.getKNearest(pos, 1);

	if (closest.size() == 0)
		return null;

	final double[] a = new double[3];
	final double[] b = new double[3];
	pos.localize(a);
	closest.get(0).getPosition().localize(b);
	final double distance = LinAlgHelpers.distance(a, b);

	if (distance <= maxDistance)
		return closest.get(0);

	return null;
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
protected void paint( final int x1, final int y1, final int x2, final int y2 )
{
	final double[] p1 = { x1, y1 };
	final double[] p2 = { x2, y2 };
	LinAlgHelpers.subtract( p2, p1, p2 );

	final double l = LinAlgHelpers.length( p2 );
	LinAlgHelpers.normalize( p2 );

	System.out.println( x1 + " " + y1 + ", " + x2 + " " + y2 + ", " + l );
	long xOld = Math.round( p1[ 0 ] ), yOld = Math.round( p1[ 1 ] );
	for ( int i = 1; i < l; ++i )
	{

		LinAlgHelpers.add( p1, p2, p1 );
		final long x = Math.round( p1[ 0 ] ), y = Math.round( p1[ 1 ] );
		if ( x != xOld || y != yOld )
		{
			paint( ( int ) x, ( int ) y );
			xOld = x;
			yOld = y;
		}
	}
	paint( x2, y2 );
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
protected void paint( final int x1, final int y1, final int x2, final int y2 )
{
	setCoordinates( x1, y1 );
	final double[] p1 = new double[ 3 ];
	final RealPoint rp1 = RealPoint.wrap( p1 );
	labelLocation.localize( p1 );

	setCoordinates( x2, y2 );
	final double[] d = new double[ 3 ];
	labelLocation.localize( d );

	LinAlgHelpers.subtract( d, p1, d );

	final double l = LinAlgHelpers.length( d );
	LinAlgHelpers.normalize( d );

	for ( int i = 1; i < l; ++i )
	{
		LinAlgHelpers.add( p1, d, p1 );
		paint( rp1 );
	}
	paint( labelLocation );
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
/**
 * Returns a boolean array indicating whether 
 * Does not modify the candidates array.
 * @param e
 * @param candidates
 * @return 
 */
public boolean[] similarlyOriented(Edgel e, List<Edgel> candidates)
{
	boolean[] match = new boolean[ candidates.size() ];
	Arrays.fill(match, false);
	
	for(int i=candidates.size()-1; i>=0; i--)
	{
		double dot = LinAlgHelpers.dot(e.getGradient(), candidates.get(i).getGradient());
		logger.trace(" dot " + dot);
		if(  dot >= 0 ){
			match[i] = true;
		}
	}
	return match;
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
/**
 * Computes an affinity between two edgels using only geometric information.
 * <P>
 * dist(e1,e1) + lam * dot
 * @param e1 first edgel
 * @param e2 second edgel
 * @param lam contribution of normal vector
 * @return
 */

public double edgelAffinitiesGeom ( Edgel e1, Edgel e2, double lam )
{
	// negate because normals should point in opposite directions
	double dot = -LinAlgHelpers.dot( e1.getGradient(), e2.getGradient() );
	
	double[] p1 = new double[e1.numDimensions()];
	e1.localize(p1);
	double[] p2 = new double[e1.numDimensions()];
	e2.localize(p2);
	
	double aff = LinAlgHelpers.squareLength( 
			ArrayUtil.subtract( p1, p2 ) );
	
	return ( 1 - lam ) * aff + ( lam * aff ) * dot ;
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
@Override
public AffineTransform3D get( final double t )
{
	final double[] qDiffCurrent = new double[ 4 ];
	final double[] qCurrent = new double[ 4 ];
	LinAlgHelpers.quaternionPower( qDiff, t, qDiffCurrent );
	LinAlgHelpers.quaternionMultiply( qStart, qDiffCurrent, qCurrent );

	final double scaleCurrent = scaleStart + t * scaleDiff;

	final double[] xg0Current = new double[ 3 ];
	final double[] tCurrent = new double[ 3 ];
	LinAlgHelpers.scale( xg0Diff, -( t * scaleEnd / scaleCurrent ), xg0Current );
	for ( int r = 0; r < 3; ++r )
		xg0Current[ r ] -= xg0Start[ r ];
	final double[][] Rcurrent = new double[ 3 ][ 3 ];
	LinAlgHelpers.quaternionToR( qCurrent, Rcurrent );
	LinAlgHelpers.mult( Rcurrent, xg0Current, tCurrent );

	final double[][] m = new double[ 3 ][ 4 ];
	for ( int r = 0; r < 3; ++r )
	{
		for ( int c = 0; c < 3; ++c )
			m[ r ][ c ] = scaleCurrent * Rcurrent[ r ][ c ];
		m[ r ][ 3 ] = scaleCurrent * tCurrent[ r ];
	}

	final AffineTransform3D transform = new AffineTransform3D();
	transform.set( m );
	return transform;
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
@Override
public AffineTransform3D get( final double t )
{

	final double scaleCurrent = scaleStart + t * scaleDiff;
	final double angleCurrent = startAngle + t * totalAngle;
	
	final double[] xg0Current = new double[ 3 ];
	final double[] tCurrent = new double[ 3 ];
	LinAlgHelpers.scale( xg0Diff, -( t * scaleEnd / scaleCurrent ), xg0Current );
	for ( int r = 0; r < 3; ++r )
		xg0Current[ r ] -= xg0Start[ r ];
	
	double[][] Rcurrent = Rotation2DHelpers.rotationToMatrix( angleCurrent );
	LinAlgHelpers.mult( Rcurrent, xg0Current, tCurrent );
	

	final double[][] m = new double[ 3 ][ 4 ];
	for ( int r = 0; r < 3; ++r )
	{
		for ( int c = 0; c < 3; ++c )
			m[ r ][ c ] = scaleCurrent * Rcurrent[ r ][ c ];
		m[ r ][ 3 ] = scaleCurrent * tCurrent[ r ];
	}
	
	final AffineTransform3D transform = new AffineTransform3D();
	transform.set( m );
	transform.set( 1.0, 2, 2 );
	
	return transform;
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
/**
 * Removes an edgel from the candidate list if its normal direction 
 * points in the opposite direction as the reference Edgel e. 
 * ( i.e. if scalar product with ref is negative ).  
 * @param e
 * @param candidates
 */
public void filterEdgelsByNormal(Edgel e, List<Edgel> candidates)
{
	
	for(int i=candidates.size()-1; i>=0; i--)
	{
		double dot = LinAlgHelpers.dot(e.getGradient(), candidates.get(i).getGradient());
		logger.trace(" dot " + dot);
		if(  dot >= 0 )
			candidates.remove(i);
	}
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
/**
	 * Returns a the dot product of the orientation of edgel e
	 * and each of the edgels in the input list
	 * @param  or orientation
	 * @param candidates
	 * @return 
	 */
	public double[] dotProduct( double[] or, Collection<Edgel> candidates)
	{
		double[] dots = new double[ candidates.size() ];
		
//		for(int i=candidates.size()-1; i>=0; i--)
		int i = 0;
		for( Edgel c : candidates )
		{
			double dot = LinAlgHelpers.dot( or, c.getGradient());
			dots[i] = dot;
			i++;
		}
		return dots;
	}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
double alignStep(final RandomAccessibleInterval< T > image, ExecutorService service)
{
	// compute error image = warped image - template
	computeDifference( Views.extendBorder( image ), currentTransform, template, error, service, Runtime.getRuntime().availableProcessors() * 2 );

	// compute transform parameter update
	final double[] gradient = new double[numParameters];

	ArrayList< Callable< Void > > calls =  new ArrayList< Callable<Void> >();
	for ( int p = 0; p < numParameters; ++p )
	{
		final int pInner = p;
		Callable< Void > callable = new Callable< Void >()
		{
			@Override
			public Void call() throws Exception
			{
				double gradT = 0;
				final Cursor< FloatType > err = Views.flatIterable( error ).cursor();
				for ( final FloatType t : Views.flatIterable( Views.hyperSlice( descent, n, pInner ) ) )
					gradT += t.getRealDouble() * err.next().getRealDouble();
				gradient[pInner] = gradT;
				return null;

			}
		};

		calls.add( callable );

	}

	List<Future<Void>> futures = null;

	try
	{
		futures = service.invokeAll( calls );
		for (Future<Void> f : futures)
			f.get();
	}
	catch ( InterruptedException | ExecutionException e)
	{
		e.printStackTrace();
	}

	final double[] dp = new double[numParameters];
	LinAlgHelpers.mult( Hinv, gradient, dp );

	// udpate transform
	currentTransform.preConcatenate( warpFunction.getAffine( dp ) );

	// return norm of parameter update vector
	return LinAlgHelpers.length( dp );
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
public synchronized void rotateView2d( boolean isClockwise )
	{
		if ( !transformEnabled )
			return;

		final SourceState< ? > source = state.getSources().get( state.getCurrentSource() );
		final AffineTransform3D sourceTransform = new AffineTransform3D();
		source.getSpimSource().getSourceTransform( state.getCurrentTimepoint(), 0, sourceTransform );

		final AffineTransform3D transform = display.getTransformEventHandler().getTransform();

		// avoid computing angle explicitly ( and thus avoid expensive inverse tan op )
		// and instead do a verbose, but faster if statements
		double m00 = transform.get( 0, 0 );
		double m01 = transform.get( 0, 1 );
		double m10 = transform.get( 1, 0 );
		double m11 = transform.get( 1, 1 );

		boolean xpos = ( m00 + m10 > 0 );
		boolean ypos = ( m01 + m11 > 0 );
		
		RotatePlane2d qTarget;
		if( isClockwise )
		{
			if( xpos && ypos )
				qTarget = RotatePlane2d.qpY;
			else if( xpos && !ypos )
				qTarget = RotatePlane2d.qnX;
			else if( !xpos && ypos )
				qTarget = RotatePlane2d.qpX;
			else if( !xpos && !ypos )
				qTarget = RotatePlane2d.qnY;
			else
				qTarget = null; 
		}
		else
		{
			if( xpos && ypos )
				qTarget = RotatePlane2d.qnY;
			else if( xpos && !ypos )
				qTarget = RotatePlane2d.qpX;
			else if( !xpos && ypos )
				qTarget = RotatePlane2d.qnX;
			else if( !xpos && !ypos )
				qTarget = RotatePlane2d.qpY;
			else
				qTarget = null; 
		}
		
		if( qTarget == null ) return;
		
		double[][] R = new double[4][4];
		LinAlgHelpers.quaternionToR( qTarget.qAlign, R );
		R[3][3] = 1.0;
		destXfm.set( R );
		
		double centerX;
		double centerY;
		if ( mouseCoordinates.isMouseInsidePanel() )
		{
			centerX = mouseCoordinates.getX();
			centerY = mouseCoordinates.getY();
		}
		else
		{
			centerY = getHeight() / 2.0;
			centerX = getWidth() / 2.0;
		}
		
//		currentAnimator = new RotationAnimator( transform, centerX, centerY, q, 300 );
//		currentAnimator = new RotationAnimator( transform, centerX, centerY, qTarget.qAlign, 300 );
		currentAnimator = new RotationAnimator2D( transform, centerX, centerY, destXfm, 300 );

		currentAnimator.setTime( System.currentTimeMillis() );
		transformChanged( transform );
	}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
public SimilarityTransformAnimator3D( final AffineTransform3D transformStart, final AffineTransform3D transformEnd, final double cX, final double cY, final long duration )
{
	super( duration );
	this.cX = cX;
	this.cY = cY;

	qStart = new double[ 4 ];
	final double[] qStartInv = new double[ 4 ];
	final double[] qEnd = new double[ 4 ];
	final double[] qEndInv = new double[ 4 ];
	qDiff = new double[ 4 ];
	Affine3DHelpers.extractRotation( transformStart, qStart );
	LinAlgHelpers.quaternionInvert( qStart, qStartInv );

	Affine3DHelpers.extractRotation( transformEnd, qEnd );
	LinAlgHelpers.quaternionInvert( qEnd, qEndInv );

	LinAlgHelpers.quaternionMultiply( qStartInv, qEnd, qDiff );
	if ( qDiff[ 0 ] < 0 )
		LinAlgHelpers.scale( qDiff, -1, qDiff );

	scaleStart = Affine3DHelpers.extractScale( transformStart, 0 );
	scaleEnd = Affine3DHelpers.extractScale( transformEnd, 0 );
	scaleDiff = scaleEnd - scaleStart;

	final double[] tStart = new double[ 3 ];
	final double[] tEnd = new double[ 3 ];
	for ( int d = 0; d < 3; ++d )
	{
		tStart[ d ] = transformStart.get( d, 3 ) / scaleStart;
		tEnd[ d ] = transformEnd.get( d, 3 ) / scaleEnd;
	}

	xg0Start = new double[3];
	final double[] xg0End = new double[3];
	xg0Diff = new double[3];

	final double[][] R = new double[ 3 ][ 3 ];
	LinAlgHelpers.quaternionToR( qStartInv, R );
	LinAlgHelpers.mult( R, tStart, xg0Start );
	LinAlgHelpers.scale( xg0Start, -1, xg0Start );
	LinAlgHelpers.quaternionToR( qEndInv, R );
	LinAlgHelpers.mult( R, tEnd, xg0End );
	LinAlgHelpers.scale( xg0End, -1, xg0End );
	LinAlgHelpers.subtract( xg0End, xg0Start, xg0Diff );
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
public SimilarityTransformAnimator2D( final AffineTransform3D transformStart, final AffineTransform3D transformEnd, final double cX, final double cY, final long duration )
{
	super( duration );
	
	startAngle = Rotation2DHelpers.extractRotation( transformStart );
	final double endAngle = Rotation2DHelpers.extractRotation( transformEnd );
	
	scaleStart = Rotation2DHelpers.extractScale( transformStart, 0 );
	scaleEnd   = Rotation2DHelpers.extractScale( transformEnd, 0 );
	scaleDiff = scaleEnd - scaleStart;
	
	totalAngle = Rotation2DHelpers.shorterAngleBetweenRotations( startAngle, endAngle );
	
	cosInvAngleStart = Math.cos( -startAngle );
	sinInvAngleStart = Math.sin( -startAngle );
	
	cosInvAngleEnd = Math.cos( -endAngle );
	sinInvAngleEnd = Math.sin( -endAngle );
	
	final double[] tStart = new double[ 3 ];
	final double[] tEnd = new double[ 3 ];
	for ( int d = 0; d < 3; ++d )
	{
		tStart[ d ] = transformStart.get( d, 3 ) / scaleStart;
		tEnd[ d ] = transformEnd.get( d, 3 ) / scaleEnd;
	}
	
	xg0Start = new double[3];
	final double[] xg0End = new double[3];
	xg0Diff = new double[3];
	
	// inverse rotate scaled starting 
	Rotation2DHelpers.cosSinMult( tStart, xg0Start, cosInvAngleStart, sinInvAngleStart );
	LinAlgHelpers.scale( xg0Start, -1, xg0Start );
		
	Rotation2DHelpers.cosSinMult( tEnd, xg0End, cosInvAngleEnd, sinInvAngleEnd );
	LinAlgHelpers.scale( xg0End, -1, xg0End );
	
	LinAlgHelpers.subtract( xg0End, xg0Start, xg0Diff );
	
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
public static void main( final String[] args )
{
	final int w = 400;
	final int h = 400;
	final int nPoints = 100000;
	final Random rand = new Random( 123124 );

	final ArrayList< Point > points = new ArrayList< Point >();
	for ( int i = 0; i < nPoints; ++i )
	{
		final long x = rand.nextInt( w );
		final long y = rand.nextInt( h );
		points.add( new Point( x, y ) );
	}

	final double[][] planes = new double[ 5 ][ 3 ]; // unit normal x, y; d

	double[] plane = planes[ 0 ];
	plane[ 0 ] = 1;
	plane[ 1 ] = 1;
	LinAlgHelpers.scale( plane, 1.0 / LinAlgHelpers.length( plane ), plane );
	plane[ 2 ] = 230;

	plane = planes[ 1 ];
	plane[ 0 ] = -1;
	plane[ 1 ] = 1;
	LinAlgHelpers.scale( plane, 1.0 / LinAlgHelpers.length( plane ), plane );
	plane[ 2 ] = -30;

	plane = planes[ 2 ];
	plane[ 0 ] = 0.1;
	plane[ 1 ] = -1;
	LinAlgHelpers.scale( plane, 1.0 / LinAlgHelpers.length( plane ), plane );
	plane[ 2 ] = -230;

	plane = planes[ 3 ];
	plane[ 0 ] = -0.5;
	plane[ 1 ] = -1;
	LinAlgHelpers.scale( plane, 1.0 / LinAlgHelpers.length( plane ), plane );
	plane[ 2 ] = -290;

	plane = planes[ 4 ];
	plane[ 0 ] = -1;
	plane[ 1 ] = 0.1;
	LinAlgHelpers.scale( plane, 1.0 / LinAlgHelpers.length( plane ), plane );
	plane[ 2 ] = -200;

	System.out.println( "partitioning list of points:" );
	BenchmarkHelper.benchmarkAndPrint( 20, false, new Runnable()
	{
		@Override
		public void run()
		{
			for ( int i = 0; i < 500; ++i )
			{
				final ArrayList< Point >[] insideoutside = getInsidePoints( points, planes );
				if ( insideoutside[ 0 ].size() > 1000000 )
					System.out.println( "bla" );
			}
		}
	} );

	System.out.println( "partitioning kdtree of points:" );
	final KDTree< Point > kdtree = new KDTree< Point >( points, points );
	final ClipConvexPolytopeKDTree< Point > clipper = new ClipConvexPolytopeKDTree< Point >( kdtree );
	BenchmarkHelper.benchmarkAndPrint( 20, false, new Runnable()
	{
		@Override
		public void run()
		{
			for ( int i = 0; i < 500; ++i )
			{
				clipper.clip( planes );
			}
		}
	} );
}
 

import net.imglib2.util.LinAlgHelpers; //导入依赖的package包/类
public static void ball() throws ImgLibException{
		String edgelMaskPrefix = "/groups/jain/home/bogovicj/projects/crackSegmentation/edgelValidate/edgelball_1";
		String edgelPrefix = "/groups/jain/home/bogovicj/projects/crackSegmentation/edgelValidate/ball_1";
		
		ImagePlusImg<FloatType,?> img = ImagePlusImgs.floats(new long[]{200, 200, 200});
		Cursor<FloatType> c = Views.iterable(Views.offset(img, new long[]{100, 100, 100})).localizingCursor();
		double[] x = new double[]{0, 0, 0};
		while (c.hasNext()) {
		  FloatType t = c.next();
		  c.localize(x);
		  if (LinAlgHelpers.squareLength(x) < 2500)
		         t.set(0xff);
		}
//		img.getImagePlus().show();
		
		int[] patchSize = new int[]{7,7,3};
		CrackCorrection<FloatType> cc = new CrackCorrection<FloatType>(
				img, img, patchSize);
		
		cc.computeEdgels();
		
		int N = 40;
		
		ArrayList<Edgel> edgels = cc.getEdgels();
		Collections.shuffle(edgels, new Random(1));
		
		ArrayImgFactory<UnsignedByteType> ubfactory = new ArrayImgFactory<UnsignedByteType>();
		Img<UnsignedByteType> edgelmaskimg = ubfactory.create(img, new UnsignedByteType(0));
		
		int i = 0;
		
		int[] patchMidPt = PatchTools.patchSizeToMidpt(patchSize);
		
		while(i < N){
			
			
			Edgel edgel = edgels.get(i);
			
//			IntType val = new IntType( i + 1 );
			UnsignedByteType val = new UnsignedByteType( 255 - i  );
			
//			AffineTransform3D xfmIn = CrackCorrection.pickTransformation(edgel);
//			CrackCorrection.edgelToView(edgel, edgelmaskimg, patchSize);
			
			
			AffineTransform3D xfmIn = EdgelTools.edgelToXfm(edgel, patchMidPt);
			double[] pos = new double[edgel.numDimensions() ];
			edgel.localize(pos);
			CrackCorrection.setMask( pos, patchSize, xfmIn, 
					Views.extendValue(edgelmaskimg, new UnsignedByteType()), val);

			i++;
		}
		
		
		ImgOps.write(edgelmaskimg, edgelMaskPrefix + ".tif");
		ImgOps.write(img, edgelPrefix + ".tif");
 
		
//		ImgUtil.write(mask, crackMaskRewriteFn);
		
	}