Skip to content

Commit b3b9799

Browse files
mfedderlymfedderly
committed
Simplest turf-isolines rewrite, update tests to check out the geojson diffs
1 parent f7f23f9 commit b3b9799

File tree

5 files changed

+857
-802
lines changed

5 files changed

+857
-802
lines changed

‎packages/turf-isolines/index.ts‎

Lines changed: 213 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -108,17 +108,226 @@ function createIsoLines(
108108

109109
const properties = { ...commonProperties, ...breaksProperties[i] };
110110
properties[zProperty] = threshold;
111-
// Pass options to marchingsquares lib to reproduce historical turf
112-
// behaviour.
113111
const isoline = multiLineString(isoContours(matrix, threshold), properties);
114112

115113
results.push(isoline);
116114
}
117115
return results;
118116
}
119117

120-
function isoContours(_matrix: number[][], _threshold: number): Position[][] {
121-
return [];
118+
function isoContours(matrix: number[][], threshold: number): Position[][] {
119+
const segments: [Position, Position][] = [];
120+
121+
for (let y = 0; y < matrix.length - 1; y++) {
122+
const rowLen = matrix[y].length - 1;
123+
for (let x = 0; x < rowLen; x++) {
124+
const tr = matrix[y][x + 1];
125+
const br = matrix[y + 1][x + 1];
126+
const bl = matrix[y + 1][x];
127+
const tl = matrix[y][x];
128+
129+
const grid =
130+
(tl >= threshold ? 8 : 0) |
131+
(tr >= threshold ? 4 : 0) |
132+
(br >= threshold ? 2 : 0) |
133+
(bl >= threshold ? 1 : 0);
134+
135+
// see https://en.wikipedia.org/wiki/Marching_squares
136+
// the two points in the line segments need to be carefully selected to keep things clockwise
137+
// and allow for efficient contour assembly
138+
switch (grid) {
139+
case 0:
140+
continue;
141+
case 1:
142+
segments.push([
143+
[x, y + frac(tl, bl)],
144+
[x + frac(bl, br), y + 1],
145+
]);
146+
break;
147+
case 2:
148+
segments.push([
149+
[x + frac(bl, br), y + 1],
150+
[x + 1, y + frac(tr, br)],
151+
]);
152+
break;
153+
case 3:
154+
segments.push([
155+
[x, y + frac(tl, bl)],
156+
[x + 1, y + frac(tr, br)],
157+
]);
158+
break;
159+
case 4:
160+
segments.push([
161+
[x + 1, y + frac(tr, br)],
162+
[x + frac(tl, tr), y],
163+
]);
164+
break;
165+
case 5: {
166+
// use the average of the 4 corners to differentiate the saddle case and correctly honor the clockwise winding
167+
const avg = (tl + tr + br + bl) / 4;
168+
const above = avg >= threshold;
169+
170+
// TODO We could divide our cell into 4 sub-cells to interpolate a 3rd middle point for each of these segments
171+
172+
if (above) {
173+
segments.push(
174+
[
175+
[x, y + frac(tl, bl)],
176+
[x + frac(tl, tr), y],
177+
],
178+
[
179+
[x + 1, y + frac(tr, br)],
180+
[x + frac(bl, br), y + 1],
181+
]
182+
);
183+
} else {
184+
segments.push(
185+
[
186+
[x + 1, y + frac(tr, br)],
187+
[x + frac(tl, tr), y],
188+
],
189+
[
190+
[x, y + frac(tl, bl)],
191+
[x + frac(bl, br), y + 1],
192+
]
193+
);
194+
}
195+
break;
196+
}
197+
case 6:
198+
segments.push([
199+
[x + frac(bl, br), y + 1],
200+
[x + frac(tl, tr), y],
201+
]);
202+
break;
203+
case 7:
204+
segments.push([
205+
[x, y + frac(tl, bl)],
206+
[x + frac(tl, tr), y],
207+
]);
208+
break;
209+
case 8:
210+
segments.push([
211+
[x + frac(tl, tr), y],
212+
[x, y + frac(tl, bl)],
213+
]);
214+
break;
215+
case 9:
216+
segments.push([
217+
[x + frac(tl, tr), y],
218+
[x + frac(bl, br), y + 1],
219+
]);
220+
break;
221+
case 10: {
222+
const avg = (tl + tr + br + bl) / 4;
223+
const above = avg >= threshold;
224+
225+
if (above) {
226+
segments.push(
227+
[
228+
[x + frac(tl, tr), y],
229+
[x + 1, y + frac(tr, br)],
230+
],
231+
[
232+
[x + frac(bl, br), y + 1],
233+
[x, y + frac(tl, bl)],
234+
]
235+
);
236+
} else {
237+
segments.push(
238+
[
239+
[x + frac(tl, tr), y],
240+
[x, y + frac(tl, bl)],
241+
],
242+
[
243+
[x + frac(bl, br), y + 1],
244+
[x + 1, y + frac(tr, br)],
245+
]
246+
);
247+
}
248+
break;
249+
}
250+
case 11:
251+
segments.push([
252+
[x + frac(tl, tr), y],
253+
[x + 1, y + frac(tr, br)],
254+
]);
255+
break;
256+
case 12:
257+
segments.push([
258+
[x + 1, y + frac(tr, br)],
259+
[x, y + frac(tl, bl)],
260+
]);
261+
break;
262+
case 13:
263+
segments.push([
264+
[x + 1, y + frac(tr, br)],
265+
[x + frac(bl, br), y + 1],
266+
]);
267+
break;
268+
case 14:
269+
segments.push([
270+
[x + frac(bl, br), y + 1],
271+
[x, y + frac(tl, bl)],
272+
]);
273+
break;
274+
case 15:
275+
// all above
276+
continue;
277+
}
278+
}
279+
}
280+
281+
// TODO we can keep track of 'active' contours while gridding above and assign line segments to a contour as we go
282+
// but then we'd need another step where we combine contours that might have been split by the saddle cases
283+
284+
const contours: Position[][] = [];
285+
286+
while (segments.length > 0) {
287+
const contour: Position[] = [...segments.shift()!];
288+
contours.push(contour);
289+
290+
let found: boolean;
291+
do {
292+
found = false;
293+
for (let i = 0; i < segments.length; i++) {
294+
const segment = segments[i];
295+
// add to the end of our contour
296+
if (
297+
segment[0][0] === contour[contour.length - 1][0] &&
298+
segment[0][1] === contour[contour.length - 1][1]
299+
) {
300+
found = true;
301+
contour.push(segment[1]);
302+
segments.splice(i, 1);
303+
break;
304+
}
305+
// add to the beginning of our contour
306+
if (
307+
segment[1][0] === contour[0][0] &&
308+
segment[1][1] === contour[0][1]
309+
) {
310+
found = true;
311+
contour.unshift(segment[0]);
312+
segments.splice(i, 1);
313+
break;
314+
}
315+
}
316+
} while (found);
317+
}
318+
319+
return contours;
320+
321+
// get the linear interpolation fraction of how far z is between z0 and z1
322+
// See https://github.com/fschutt/marching-squares/blob/master/src/lib.rs
323+
function frac(z0: number, z1: number): number {
324+
if (z0 === z1) {
325+
return 0.5;
326+
}
327+
328+
let t = (threshold - z0) / (z1 - z0);
329+
return t > 1 ? 1 : t < 0 ? 0 : t;
330+
}
122331
}
123332

124333
/**

0 commit comments

Comments
 (0)