/
GaussianLineSmoothing.java
144 lines (128 loc) · 4.62 KB
/
GaussianLineSmoothing.java
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
package eu.europa.ec.eurostat.jgiscotools.algo.line;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.LineString;
import eu.europa.ec.eurostat.jgiscotools.algo.base.DouglasPeuckerRamerFilter;
/**
* Line gaussian smoothing.
*
* @author julien Gaffuri
*
*/
public class GaussianLineSmoothing {
/**
* @param line
* @param sigmaM
* @return
*/
public static LineString get(LineString line, double sigmaM){ return get(line, sigmaM, -1); }
/**
* Line gaussian smoothing.
* The position of each point is the average position of its neighbors, weighted by a gaussian kernel.
* For non-closed lines, the initial and final points are preserved.
*
* @param line The input line
* @param sigmaM The standard deviation of the gaussian kernel. The larger, the more smoothed.
* @param resolution The target resolution of the geometry. This parameter is used to filter the final geometry.
* @return
*/
public static LineString get(LineString line, double sigmaM, double resolution){
if(line.getCoordinates().length <= 2) return (LineString) line.copy();
boolean isClosed = line.isClosed();
double length = line.getLength();
double densifiedResolution = sigmaM/3;
//handle extreme cases: too large sigma resulting in too large densified resolution.
if(densifiedResolution > 0.25*length ) {
if(isClosed){
//return tiny triangle nearby center point
Coordinate c = line.getCentroid().getCoordinate();
length *= 0.01;
return line.getFactory().createLineString(new Coordinate[]{ new Coordinate(c.x-length,c.y-length), new Coordinate(c.x,c.y+length), new Coordinate(c.x+length,c.y-length), new Coordinate(c.x-length,c.y-length) });
} else {
//return segment
return line.getFactory().createLineString(new Coordinate[]{ line.getCoordinateN(0), line.getCoordinateN(line.getNumPoints()-1) });
}
}
//compute densified line
Coordinate[] densifiedCoords = LittleThumblingDensifier.densify(line, densifiedResolution).getCoordinates();
//build ouput line structure
int nb = (int) (length/densifiedResolution);
Coordinate[] out = new Coordinate[nb+1];
//prepare gaussian coefficients
int n = 7*3; //it should be: E(7*sigma/densifiedResolution), which is 7*3;
double gcs[] = new double[n+1];
{
double a = sigmaM*Math.sqrt(2*Math.PI);
double b = sigmaM*sigmaM*2;
double d = densifiedResolution*densifiedResolution;
for(int i=0; i<n+1; i++) gcs[i] = Math.exp(-i*i*d/b) /a;
}
Coordinate c0 = densifiedCoords[0];
Coordinate cN = densifiedCoords[nb];
for(int i=0; i<nb; i++) {
if(!isClosed && i==0) continue;
//compute coordinates of point i of the smoothed line (gauss mean)
double x=0.0, y=0.0;
for(int j=-n; j<=n; j++) {
//index of the point to consider on the original densified line
int q = i+j;
//find coordinates (xq,yq) of point q
double xq, yq;
if(q<0) {
if(isClosed) {
//make loop to get the right point
q = q%nb; if(q<0) q+=nb;
Coordinate c = densifiedCoords[q];
xq = c.x;
yq = c.y;
} else {
//get symetric point
q = (-q)%nb; if(q==0) q=nb;
Coordinate c = densifiedCoords[q];
xq = 2*c0.x-c.x;
yq = 2*c0.y-c.y;
}
} else if (q>nb) {
if(isClosed) {
//make loop to get the right point
q = q%nb; if(q==0) q=nb;
Coordinate c = densifiedCoords[q];
xq = c.x;
yq = c.y;
} else {
//get symetric point
q = nb-q%nb; if(q==nb) q=0;
Coordinate c = densifiedCoords[q];
xq = 2*cN.x-c.x;
yq = 2*cN.y-c.y;
}
} else {
//general case (most frequent)
Coordinate c = densifiedCoords[q];
xq = c.x;
yq = c.y;
}
//get gaussian coefficient
double gc = gcs[j>=0?j:-j];
//add contribution of point q to new position of point i
x += xq*gc;
y += yq*gc;
}
//assign smoothed position of point i
out[i] = new Coordinate(x*densifiedResolution, y*densifiedResolution);
}
//handle start and end points
if(isClosed) {
//ensure start and end locations are the same
out[nb] = out[0];
} else {
//ensure start and end points are at the same position as the initial geometry
out[0] = densifiedCoords[0];
out[nb] = densifiedCoords[densifiedCoords.length-1];
}
//prepare final line, applying some filtering
LineString lsOut = line.getFactory().createLineString(out);
if(resolution<0) resolution = densifiedResolution /3;
lsOut = (LineString) DouglasPeuckerRamerFilter.get( lsOut , resolution);
return lsOut;
}
}