use of org.opengis.geometry.MismatchedDimensionException in project sis by apache.
the class Line method fit.
/**
* Given a sequence of points, fits them to a straight line
* <var>y</var> = <var>slope</var>⋅<var>x</var> + <var>y₀</var> in a least-squares senses.
* Points shall be two dimensional with ordinate values in the (<var>x</var>,<var>y</var>) order.
* This method assumes that the <var>x</var> values are precise and all uncertainty is in <var>y</var>.
* {@link Double#NaN} ordinate values are ignored.
*
* @param points the two-dimensional points.
* @return estimation of the correlation coefficient. The closer this coefficient is to +1 or -1, the better the fit.
* @throws MismatchedDimensionException if a point is not two-dimensional.
*/
public double fit(final Iterable<? extends DirectPosition> points) {
ArgumentChecks.ensureNonNull("points", points);
int i = 0, n = 0;
final DoubleDouble mean_x = new DoubleDouble();
final DoubleDouble mean_y = new DoubleDouble();
for (final DirectPosition p : points) {
final int dimension = p.getDimension();
if (dimension != DIMENSION) {
throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_3, "points[" + i + ']', DIMENSION, dimension));
}
i++;
final double x, y;
if (// Test first the dimension which is most likely to contain NaN.
!isNaN(y = p.getOrdinate(1)) && !isNaN(x = p.getOrdinate(0))) {
mean_x.add(x);
mean_y.add(y);
n++;
}
}
mean_x.divide(n, 0);
mean_y.divide(n, 0);
/*
* We have to solve two equations with two unknowns:
*
* 1) mean(y) = m⋅mean(x) + y₀
* 2) mean(x⋅y) = m⋅mean(x²) + y₀⋅mean(x)
*
* Those formulas lead to a quadratic equation. However, the formulas become very simples
* if we set 'mean(x) = 0'. We can achieve this result by computing instead of (2):
*
* 2b) mean(Δx⋅y) = m⋅mean(Δx²)
*
* where dx = x-mean(x). In this case mean(Δx) == 0.
*/
final DoubleDouble a = new DoubleDouble();
final DoubleDouble b = new DoubleDouble();
final DoubleDouble mean_x2 = new DoubleDouble();
final DoubleDouble mean_y2 = new DoubleDouble();
final DoubleDouble mean_xy = new DoubleDouble();
for (final DirectPosition p : points) {
final double y;
if (// Test first the dimension which is most likely to contain NaN.
!isNaN(y = p.getOrdinate(1)) && !isNaN(a.value = p.getOrdinate(0))) {
a.error = 0;
// Δx = x - mean_x
a.subtract(mean_x);
b.setFrom(a);
// (Δx)² = Δx * Δx
b.multiply(a);
// mean_x² += (Δx)²
mean_x2.add(b);
// xy = Δx * y
a.multiply(y);
// mean_xy += xy
mean_xy.add(a);
// y² = y * y
a.setToProduct(y, y);
// mean_y² += y²
mean_y2.add(a);
}
}
mean_x2.divide(n, 0);
mean_y2.divide(n, 0);
mean_xy.divide(n, 0);
/*
* Assuming that 'mean(x) == 0', then the correlation
* coefficient can be approximate by:
*
* R = mean(xy) / sqrt( mean(x²) * (mean(y²) - mean(y)²) )
*/
a.setFrom(mean_xy);
// slope = mean_xy / mean_x²
a.divide(mean_x2);
b.setFrom(mean_x);
b.multiply(a);
b.negate();
// y₀ = mean_y - mean_x * slope
b.add(mean_y);
setEquation(a, b);
/*
* Compute the correlation coefficient:
* mean_xy / sqrt(mean_x2 * (mean_y2 - mean_y * mean_y))
*/
a.setFrom(mean_y);
a.multiply(mean_y);
a.negate();
a.add(mean_y2);
a.multiply(mean_x2);
a.sqrt();
a.inverseDivide(mean_xy);
return a.value;
}
use of org.opengis.geometry.MismatchedDimensionException in project sis by apache.
the class ConcatenatedTransform method create.
/**
* Concatenates the two given transforms.
* If the concatenation result works with two-dimensional input and output points,
* then the returned transform will implement {@link MathTransform2D}.
* Likewise if the concatenation result works with one-dimensional input and output points,
* then the returned transform will implement {@link MathTransform1D}.
*
* <div class="note"><b>Implementation note:</b>
* {@code ConcatenatedTransform} implementations are available in two versions: direct and non-direct.
* The "non-direct" versions use an intermediate buffer when performing transformations; they are slower
* and consume more memory. They are used only as a fallback when a "direct" version can not be created.</div>
*
* @param tr1 the first math transform.
* @param tr2 the second math transform.
* @param factory the factory which is (indirectly) invoking this method, or {@code null} if none.
* @return the concatenated transform.
*
* @see MathTransforms#concatenate(MathTransform, MathTransform)
*/
public static MathTransform create(MathTransform tr1, MathTransform tr2, final MathTransformFactory factory) throws FactoryException, MismatchedDimensionException {
final int dim1 = tr1.getTargetDimensions();
final int dim2 = tr2.getSourceDimensions();
if (dim1 != dim2) {
throw new MismatchedDimensionException(Resources.format(Resources.Keys.CanNotConcatenateTransforms_2, getName(tr1), getName(tr2)) + ' ' + Errors.format(Errors.Keys.MismatchedDimension_2, dim1, dim2));
}
MathTransform mt = createOptimized(tr1, tr2, factory);
if (mt != null) {
return mt;
}
/*
* If at least one math transform is an instance of ConcatenatedTransform and assuming
* that MathTransforms are associatives, tries the following arrangements and select
* the one with the fewest amount of steps:
*
* Assuming : tr1 = (A * B)
* tr2 = (C * D)
*
* Current : (A * B) * (C * D) Will be the selected one if nothing better.
* Try k=0 : A * (B * (C * D)) Implies A * ((B * C) * D) through recursivity.
* Try k=1 : ((A * B) * C) * D Implies (A * (B * C)) * D through recursivity.
* Try k=2 : Tried only if try k=1 changed something.
*
* TODO: The same combination may be computed more than once (e.g. (B * C) above).
* Should not be a big deal if there is not two many steps. In the even where
* it would appears a performance issue, we could maintain a Map of combinations
* already computed. The map would be local to a "create" method execution.
*/
int stepCount = getStepCount(tr1) + getStepCount(tr2);
// Really 'true' because we want at least 2 iterations.
boolean tryAgain = true;
for (int k = 0; ; k++) {
MathTransform c1 = tr1;
MathTransform c2 = tr2;
final boolean first = (k & 1) == 0;
MathTransform candidate = first ? c1 : c2;
while (candidate instanceof ConcatenatedTransform) {
final ConcatenatedTransform ctr = (ConcatenatedTransform) candidate;
if (first) {
c1 = candidate = ctr.transform1;
c2 = create(ctr.transform2, c2, factory);
} else {
c1 = create(c1, ctr.transform1, factory);
c2 = candidate = ctr.transform2;
}
final int c = getStepCount(c1) + getStepCount(c2);
if (c < stepCount) {
tr1 = c1;
tr2 = c2;
stepCount = c;
tryAgain = true;
}
}
if (!tryAgain)
break;
tryAgain = false;
}
/*
* Tries again the check for optimized cases (identity, etc.), because a
* transform may have been simplified to identity as a result of the above.
*/
mt = createOptimized(tr1, tr2, factory);
if (mt != null) {
return mt;
}
/*
* Can not avoid the creation of a ConcatenatedTransform object.
* Check for the type to create (1D, 2D, general case...)
*/
final int dimSource = tr1.getSourceDimensions();
final int dimTarget = tr2.getTargetDimensions();
if (dimSource == 1 && dimTarget == 1) {
/*
* Result needs to be a MathTransform1D.
*/
if (tr1 instanceof MathTransform1D && tr2 instanceof MathTransform1D) {
return new ConcatenatedTransformDirect1D((MathTransform1D) tr1, (MathTransform1D) tr2);
} else {
return new ConcatenatedTransform1D(tr1, tr2);
}
} else if (dimSource == 2 && dimTarget == 2) {
/*
* Result needs to be a MathTransform2D.
*/
if (tr1 instanceof MathTransform2D && tr2 instanceof MathTransform2D) {
return new ConcatenatedTransformDirect2D((MathTransform2D) tr1, (MathTransform2D) tr2);
} else {
return new ConcatenatedTransform2D(tr1, tr2);
}
} else if (// dim1 = tr1.getTargetDimensions() and
dimSource == tr1.getTargetDimensions() && // dim2 = tr2.getSourceDimensions() may not be true anymore.
dimTarget == tr2.getSourceDimensions()) {
return new ConcatenatedTransformDirect(tr1, tr2);
} else {
return new ConcatenatedTransform(tr1, tr2);
}
}
use of org.opengis.geometry.MismatchedDimensionException in project mkgmap by openstreetmap.
the class PrecompSeaGenerator method runSeaGeneration.
public void runSeaGeneration() throws MismatchedDimensionException, TransformException, IOException, InterruptedException {
createShapefileAccess();
// get all tiles that need to be processed
List<uk.me.parabola.imgfmt.app.Area> remainingTiles = getTiles();
// initialize the count down so that it is possible to get the
// information when all tiles are finished
CountDownLatch tilesCountdown = new CountDownLatch(remainingTiles.size());
// start a printer that outputs how many tiles still need to be
// processed
new ProgressPrinter(tilesCountdown).start();
// start the saver thread that stores the tiles to disc and creates
// the index file
PrecompSeaSaver precompSaver = new PrecompSeaSaver(outputDir, usePbfFormat);
new Thread(precompSaver, "SaveThread").start();
// requirements
while (remainingTiles.isEmpty() == false) {
// create a list with all tiles that are processed within this cycle
List<uk.me.parabola.imgfmt.app.Area> tiles = new ArrayList<uk.me.parabola.imgfmt.app.Area>();
tiles.addAll(remainingTiles.subList(0, Math.min(tilesPerCycle, remainingTiles.size())));
remainingTiles.subList(0, Math.min(tilesPerCycle, remainingTiles.size())).clear();
// create the mergers that merge the data of one tile
List<PrecompSeaMerger> mergers = createMergers(tiles, tilesCountdown, precompSaver.getQueue());
// create an overall area for a simple check if a polygon read from the
// shape file intersects one of the currently processed sea tiles
Area tileArea = new Area();
for (PrecompSeaMerger m : mergers) {
tileArea.add(new Area(m.getTileBounds()));
// start the mergers
service.execute(m);
}
openShapefile();
int numPolygon = 0;
long lastInfo = System.currentTimeMillis();
// read all polygons from the shape file and add them to the queues of the
// merger threads
Geometry wgs84Poly = null;
while (shapeIterator.hasNext()) {
Feature feature = shapeIterator.next();
GeometryAttribute geom = feature.getDefaultGeometryProperty();
Geometry poly = (Geometry) geom.getValue();
if (poly == null) {
continue;
}
try {
wgs84Poly = transformToWGS84(poly);
} catch (Exception exp) {
System.err.println(exp);
continue;
}
if (wgs84Poly.getNumGeometries() != 1) {
// only simple polygons are supported by now
// maybe this could be changed in future?
System.err.println("Polygon from shapefile has " + wgs84Poly.getNumGeometries() + " geometries. Only one geometry is supported.");
System.err.println("Skip polygon.");
continue;
}
Geometry bounds = wgs84Poly.getEnvelope();
if (bounds.isEmpty()) {
System.err.println("Empty or non polygon: " + bounds);
} else {
Area polyBounds = convertToArea(bounds);
// currently processed
if (polyBounds.intersects(tileArea.getBounds2D())) {
// yes it touches at least one tile => convert it to
// a java.awt.geom.Area object
Area polyAsArea = convertToArea(wgs84Poly.getGeometryN(0));
// polygon to the queues of them
for (PrecompSeaMerger mThread : mergers) {
if (mThread.getTileBounds().intersects(polyAsArea.getBounds2D())) {
try {
mThread.getQueue().put(polyAsArea);
} catch (InterruptedException exp) {
exp.printStackTrace();
}
}
}
}
numPolygon++;
if ((numPolygon) % 50000 == 0 || System.currentTimeMillis() - lastInfo > 30000) {
// print out the current number of polygons already processed
System.out.println("Worked out " + (numPolygon) + " polygons");
lastInfo = System.currentTimeMillis();
}
}
}
closeShapefile();
System.out.println("Reading shapefile finished");
// signal all mergers that all polygons have been read
for (PrecompSeaMerger mThread : mergers) {
mThread.signalInputComplete();
}
// may occurr
while (tilesCountdown.getCount() > remainingTiles.size() + 2 * tilesPerCycle) {
Thread.sleep(50L);
}
}
// wait until all tiles have been merged
tilesCountdown.await();
// wait until the saver for the tiles is finished
precompSaver.waitForFinish();
// shutdown the executor service
service.shutdown();
}
Aggregations