mpicbg.models.IllDefinedDataPointsException Java Examples

/**
 * Fits sampled points to a model.
 *
 * Stolen from
 *
 * <a href="https://github.com/axtimwalde/fiji-scripts/blob/master/TrakEM2/visualize-ct-difference.bsh#L90-L106">
 *     https://github.com/axtimwalde/fiji-scripts/blob/master/TrakEM2/visualize-ct-difference.bsh#L90-L106
 * </a>.
 *
 * @param  model                model to fit (note: model will be changed by this operation).
 * @param  coordinateTransform  transform to apply to each sampled point.
 * @param  sampleWidth          width of each sample.
 * @param  sampleHeight         height of each sample.
 * @param  samplesPerDimension  number of samples to take in each dimension.
 */
public static void fit(final Model<?> model,
                       final CoordinateTransform coordinateTransform,
                       final double sampleWidth,
                       final double sampleHeight,
                       final int samplesPerDimension)
        throws NotEnoughDataPointsException, IllDefinedDataPointsException {

    final List<PointMatch> matches = new ArrayList<>();

    for (int y = 0; y < samplesPerDimension; ++y) {
        final double sampleY = y * sampleHeight;
        for (int x = 0; x < samplesPerDimension; ++x) {
            final double sampleX = x * sampleWidth;
            final Point p = new Point(new double[]{sampleX, sampleY});
            p.apply(coordinateTransform);
            matches.add(new PointMatch(p, p));
        }
    }

    model.fit(matches);
}
 

/**
 * Sample the average scaling of a given {@link CoordinateTransform} by transferring a set of point samples using
 * the {@link CoordinateTransform} and then least-squares fitting a {@link SimilarityModel2D} to it.
 *
 * @param width  of the samples set
 * @param height of the samples set
 * @param dx     spacing between samples
 *
 * @return average scale factor
 */
public static double sampleAverageScale(final CoordinateTransform ct,
                                        final int width,
                                        final int height,
                                        final double dx) {
    final ArrayList<PointMatch> samples = new ArrayList<>();
    for (double y = 0; y < height; y += dx) {
        for (double x = 0; x < width; x += dx) {
            final Point p = new Point(new double[]{x, y});
            p.apply(ct);
            samples.add(new PointMatch(p, p));
        }
    }
    final AffineModel2D model = new AffineModel2D();
    try {
        model.fit(samples);
    } catch (final NotEnoughDataPointsException | IllDefinedDataPointsException e) {
        LOG.warn("failed to fit samples, returning scale factor of 1", e);
        return 1;
    }
    final double[] data = new double[6];
    model.toArray(data);
    // return 1;
    return Math.sqrt(Math.max(data[0] * data[0] + data[1] * data[1], data[2] * data[2] + data[3] * data[3]));
}
 

/**
 * Estimates an initial {@link Model} based on some given {@link PointMatch}es. Note that the {@link PointMatch}es have to be stored as PointMatch(target,reference). 
 * 
 * @param matches - The {@link List} of apriori known {@link PointMatch}es
 * @param model - The {@link Model} to use
 * @throws NotEnoughDataPointsException
 * @throws IllDefinedDataPointsException
 */
public void estimateIntialModel( final List<PointMatch> matches, final Model<?> model ) throws NotEnoughDataPointsException, IllDefinedDataPointsException
{
	/* remove ambigous correspondences */
	ambigousMatches = removeAmbigousMatches( matches );
	
	/* fit the model */
	model.fit( matches );

	/* apply the new model of the target to determine the error */
	for ( final P point : target )
		point.apply( model );

	/* compute the output */
	avgError = PointMatch.meanDistance( matches );
	maxError = PointMatch.maxDistance( matches );
	numMatches = matches.size();
	pointMatches = matches;		
}
 

@Override
   public < P extends PointMatch >void fit(final Collection< P > matches)
		throws NotEnoughDataPointsException, IllDefinedDataPointsException
{

	final Stack< java.awt.Point > sourcePoints = new Stack<java.awt.Point>();
	final Stack< java.awt.Point > targetPoints = new Stack<java.awt.Point>();

	for ( final P pm : matches )
	{
		final double[] p1 = pm.getP1().getL();
		final double[] p2 = pm.getP2().getL();

		targetPoints.add( new java.awt.Point( ( int )Math.round( p1[ 0 ] ), ( int )Math.round( p1[ 1 ] ) ) );
		sourcePoints.add( new java.awt.Point( ( int )Math.round( p2[ 0 ] ), ( int )Math.round( p2[ 1 ] ) ) );
	}

	final Transformation transf = bUnwarpJ_.computeTransformationBatch(sourceWidth,
			sourceHeight, width, height, sourcePoints, targetPoints, parameter);
	this.set(transf.getIntervals(), transf.getDirectDeformationCoefficientsX(),
			transf.getDirectDeformationCoefficientsY(), width, height);
}
 

public void calculate(final List<List<PointMatch>> consensusSetInliers,
                      final Model aggregateModel)
        throws NotEnoughDataPointsException, IllDefinedDataPointsException {

    consensusSetDeltaXStandardDeviations = new ArrayList<>();
    consensusSetDeltaYStandardDeviations = new ArrayList<>();
    aggregateDeltaXAndYStandardDeviation = new double[] { 0.0, 0.0 };

    final List<PointMatch> aggregatedInliers = new ArrayList<>();
    for (final List<PointMatch> consensusSet : consensusSetInliers) {
        if (consensusSet.size() > 0) {
            final double[] worldDeltaXAndYStandardDeviation = getWorldDeltaXAndYStandardDeviation(consensusSet);
            consensusSetDeltaXStandardDeviations.add(worldDeltaXAndYStandardDeviation[0]);
            consensusSetDeltaYStandardDeviations.add(worldDeltaXAndYStandardDeviation[1]);
            if (aggregateModel != null) {
                consensusSet.forEach(pm -> aggregatedInliers.add(new PointMatch(pm.getP1().clone(),
                                                                                pm.getP2().clone())));
            }
        }
    }

    if (aggregateModel == null) {
        if (consensusSetDeltaXStandardDeviations.size() > 0) {
            aggregateDeltaXAndYStandardDeviation = new double[]{
                    consensusSetDeltaXStandardDeviations.get(0),
                    consensusSetDeltaYStandardDeviations.get(0)
            };
        }
    } else {
        aggregateModel.fit(aggregatedInliers);
        this.aggregateDeltaXAndYStandardDeviation = getWorldDeltaXAndYStandardDeviation(aggregatedInliers);
    }
}
 

@Override
final public <P extends PointMatch> void fit( final Collection< P > matches )
	throws NotEnoughDataPointsException, IllDefinedDataPointsException
{
	if ( matches.size() < MIN_NUM_MATCHES )
		throw new NotEnoughDataPointsException( matches.size() + " matches given, we need at least " + MIN_NUM_MATCHES + " data point." );
}
 

@Override
final public <P extends PointMatch> void fit( final Collection< P > matches )
	throws NotEnoughDataPointsException, IllDefinedDataPointsException
{
	if ( matches.size() < MIN_NUM_MATCHES )
		throw new NotEnoughDataPointsException( matches.size() + " matches given, we need at least " + MIN_NUM_MATCHES + " data point." );
}
 

@Override
public < P extends PointMatch >void fit( final Collection< P > pointMatches ) throws NotEnoughDataPointsException, IllDefinedDataPointsException
{
	if ( pointMatches.size() < getMinNumMatches() )
		throw new NotEnoughDataPointsException( pointMatches.size() + " data points are not enough to estimate a 2d polynomial of order " + nlt.getDimension() + ", at least " + getMinNumMatches() + " data points required." );

	affine.fit( pointMatches );

	final double h1[][] = new double[ pointMatches.size() ][ 2 ];
    final double h2[][] = new double[ pointMatches.size() ][ 2 ];

    int i = 0;
	for ( final P match : pointMatches )
    {
    	final double[] tmp1 = match.getP1().getL().clone();
    	affine.applyInPlace( tmp1 );

    	final double[] tmp2 = match.getP2().getW();

    	h1[ i ] = new double[]{ tmp1[ 0 ], tmp1[ 1 ] };
		h2[ i ] = new double[]{ tmp2[ 0 ], tmp2[ 1 ] };

		++i;
    }

	nlt.fit( h1, h2, lambda );
}
 

/**
 * Minimize the displacement of all {@link PointMatch Correspondence pairs}
 * of all {@link Tile Tiles}
 * 
 * @param maxAllowedError do not accept convergence if error is < max_error
 * @param maxIterations stop after that many iterations even if there was
 *   no minimum found
 * @param maxPlateauwidth convergence is reached if the average absolute
 *   slope in an interval of this size and half this size is smaller than
 *   0.0001 (in double accuracy).  This is assumed to prevent the algorithm
 *   from stopping at plateaus smaller than this value.
 * @param debugLevel defines if the Optimizer prints the output at the end of the process
 */
public void optimize(
		final double maxAllowedError,
		final int maxIterations,
		final int maxPlateauwidth,
		final int debugLevel ) throws NotEnoughDataPointsException, IllDefinedDataPointsException 
{
	final ErrorStatistic observer = new ErrorStatistic( maxPlateauwidth + 1 );
	
	int i = 0;
	
	boolean proceed = i < maxIterations;
	
	while ( proceed )
	{
		for ( final TileSPIM tile : tiles )
		{
			if ( fixedTiles.contains( tile ) ) continue;
			tile.update();
			tile.fitModel();
			tile.update();
		}
		update();
		observer.add( error );
		
		if ( i > maxPlateauwidth )
		{
			proceed = error > maxAllowedError;
			
			int d = maxPlateauwidth;
			while ( !proceed && d >= 1 )
			{
				try
				{
					proceed |= Math.abs( observer.getWideSlope( d ) ) > 0.0001;
				}
				catch ( Exception e ) { e.printStackTrace(); }
				d /= 2;
			}
		}
		
		proceed &= ++i < maxIterations;
	}
	
	if ( debugLevel <= ViewStructure.DEBUG_MAIN )
	{
		println( "Successfully optimized configuration of " + tiles.size() + " tiles after " + i + " iterations:" );
		println( "  average displacement: " + decimalFormat.format( error ) + "px" );
		println( "  minimal displacement: " + decimalFormat.format( minError ) + "px" );
		println( "  maximal displacement: " + decimalFormat.format( maxError ) + "px" );
	}
}
 

@Override
final public <P extends PointMatch> void fit( final Collection< P > matches )
	throws NotEnoughDataPointsException, IllDefinedDataPointsException
{
	if ( matches.size() < MIN_NUM_MATCHES )
		throw new NotEnoughDataPointsException( matches.size() + " data points are not enough to estimate a 3d affine model, at least " + MIN_NUM_MATCHES + " data points required." );

	double
		a00, a01, a02,
		     a11, a12,
		          a22;
	double
		b00, b01, b02,
		b10, b11, b12,
		b20, b21, b22;

	a00 = a01 = a02 = a11 = a12 = a22 = b00 = b01 = b02 = b10 = b11 = b12 = b20 = b21 = b22 = 0;
	for ( final PointMatch m : matches )
	{
		final double[] p = m.getP1().getL();
		final double[] q = m.getP2().getW();
		final double w = m.getWeight();

		final double px = p[ 0 ], py = p[ 1 ], pz = p[ 2 ];
		final double qx = q[ 0 ], qy = q[ 1 ], qz = q[ 2 ];
		a00 += w * px * px;
		a01 += w * px * py;
		a02 += w * px * pz;
		a11 += w * py * py;
		a12 += w * py * pz;
		a22 += w * pz * pz;

		b00 += w * px * qx;
		b01 += w * px * qy;
		b02 += w * px * qz;
		b10 += w * py * qx;
		b11 += w * py * qy;
		b12 += w * py * qz;
		b20 += w * pz * qx;
		b21 += w * pz * qy;
		b22 += w * pz * qz;
	}

	final double det =
		a00 * a11 * a22 +
		a01 * a12 * a02 +
		a02 * a01 * a12 -
		a02 * a11 * a02 -
		a12 * a12 * a00 -
		a22 * a01 * a01;

	if ( det == 0 )
		throw new IllDefinedDataPointsException();

	final double idet = 1f / det;

	final double ai00 = ( a11 * a22 - a12 * a12 ) * idet;
	final double ai01 = ( a02 * a12 - a01 * a22 ) * idet;
	final double ai02 = ( a01 * a12 - a02 * a11 ) * idet;
	final double ai11 = ( a00 * a22 - a02 * a02 ) * idet;
	final double ai12 = ( a02 * a01 - a00 * a12 ) * idet;
	final double ai22 = ( a00 * a11 - a01 * a01 ) * idet;

	m00 = ai00 * b00 + ai01 * b10 + ai02 * b20;
	m01 = ai01 * b00 + ai11 * b10 + ai12 * b20;
	m02 = ai02 * b00 + ai12 * b10 + ai22 * b20;

	m10 = ai00 * b01 + ai01 * b11 + ai02 * b21;
	m11 = ai01 * b01 + ai11 * b11 + ai12 * b21;
	m12 = ai02 * b01 + ai12 * b11 + ai22 * b21;

	m20 = ai00 * b02 + ai01 * b12 + ai02 * b22;
	m21 = ai01 * b02 + ai11 * b12 + ai12 * b22;
	m22 = ai02 * b02 + ai12 * b12 + ai22 * b22;
}
 

@Override
final public <P extends PointMatch> void fit( final Collection< P > matches )
	throws NotEnoughDataPointsException, IllDefinedDataPointsException
{
	if ( matches.size() < MIN_NUM_MATCHES )
		throw new NotEnoughDataPointsException( matches.size() + " data points are not enough to estimate a 3d rigid model, at least " + MIN_NUM_MATCHES + " data points required." );

	// calculate N
	double Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz;
	Sxx = Sxy = Sxz = Syx = Syy = Syz = Szx = Szy = Szz = 0;
	
	for ( final PointMatch m : matches )
	{
		final double[] p = m.getP1().getL(); 
		final double[] q = m.getP2().getW();
		
		final double x1 = p[ 0 ];
		final double y1 = p[ 1 ];
		final double z1 = p[ 2 ];
		final double x2 = q[ 0 ];
		final double y2 = q[ 1 ];
		final double z2 = q[ 2 ];
		Sxx += x1 * x2;
		Sxy += x1 * y2;
		Sxz += x1 * z2;
		Syx += y1 * x2;
		Syy += y1 * y2;
		Syz += y1 * z2;
		Szx += z1 * x2;
		Szy += z1 * y2;
		Szz += z1 * z2;
	}
	
	N[0][0] = Sxx + Syy + Szz;
	N[0][1] = Syz - Szy;
	N[0][2] = Szx - Sxz;
	N[0][3] = Sxy - Syx;
	N[1][0] = Syz - Szy;
	N[1][1] = Sxx - Syy - Szz;
	N[1][2] = Sxy + Syx;
	N[1][3] = Szx + Sxz;
	N[2][0] = Szx - Sxz;
	N[2][1] = Sxy + Syx;
	N[2][2] = -Sxx + Syy - Szz;
	N[2][3] = Syz + Szy;
	N[3][0] = Sxy - Syx;
	N[3][1] = Szx + Sxz;
	N[3][2] = Syz + Szy;
	N[3][3] = -Sxx - Syy + Szz;

	// calculate eigenvector with maximal eigenvalue
	/*
	final JacobiFloat jacobi = new JacobiFloat(N);
	final float[][] eigenvectors = jacobi.getEigenVectors();
	final float[] eigenvalues = jacobi.getEigenValues();
	int index = 0;
	for (int i = 1; i < 4; i++)
		if (eigenvalues[i] > eigenvalues[index])
			index = i;

	final float[] q = eigenvectors[index];
	final float q0 = q[0], qx = q[1], qy = q[2], qz = q[3];
	*/
	// calculate eigenvector with maximal eigenvalue

	final EigenvalueDecomposition evd = new EigenvalueDecomposition( new Matrix( N ) );
	
	final double[] eigenvalues = evd.getRealEigenvalues();
	final Matrix eigenVectors = evd.getV();

	int index = 0;
	for (int i = 1; i < 4; i++)
		if (eigenvalues[i] > eigenvalues[index])
			index = i;

	final double q0 = eigenVectors.get( 0, index ); 
	final double qx = eigenVectors.get( 1, index );
	final double qy = eigenVectors.get( 2, index );
	final double qz = eigenVectors.get( 3, index );
	
	// computational result
	
	// rotational part
	m00 = (q0 * q0 + qx * qx - qy * qy - qz * qz);
	m01 = 2 * (qx * qy - q0 * qz);
	m02 = 2 * (qx * qz + q0 * qy);
	m10 = 2 * (qy * qx + q0 * qz);
	m11 = (q0 * q0 - qx * qx + qy * qy - qz * qz);
	m12 = 2 * (qy * qz - q0 * qx);
	m20 = 2 * (qz * qx - q0 * qy);
	m21 = 2 * (qz * qy + q0 * qx);
	m22 = (q0 * q0 - qx * qx - qy * qy + qz * qz);
	
	/*
	// translational part
	result.apply(c1x, c1y, c1z);
	result.a03 = c2x - result.x;
	result.a13 = c2y - result.y;
	result.a23 = c2z - result.z;
	*/
}