conrec is a Go translation of the C version of CONREC by Paul Bourke: http://paulbourke.net/papers/conrec/conrec.c conrec takes g, an m×n grid function, a sorted slice of contour heights and a conrecLine function. For full details of the algorithm, see the paper at http://paulbourke.net/papers/con
(g GridXYZ, heights []float64, fn conrecLine)
| 33 | // For full details of the algorithm, see the paper at |
| 34 | // http://paulbourke.net/papers/conrec/ |
| 35 | func conrec(g GridXYZ, heights []float64, fn conrecLine) { |
| 36 | var ( |
| 37 | p1, p2 point |
| 38 | |
| 39 | h [5]float64 |
| 40 | sh [5]int |
| 41 | xh, yh [5]float64 |
| 42 | |
| 43 | im = [4]int{0, 1, 1, 0} |
| 44 | jm = [4]int{0, 0, 1, 1} |
| 45 | |
| 46 | // We differ from conrec.c in the assignment of a single value |
| 47 | // in cases (castab in conrec.c). The value of castab[1][1][1] is |
| 48 | // 3, but we set cases[1][1][1] to 0. |
| 49 | // |
| 50 | // axiom: When we have a section of the grid where all the |
| 51 | // Z values are equal, and equal to a contour height we would |
| 52 | // expect to have no internal segments to draw. |
| 53 | // |
| 54 | // This is covered by case g) in Paul Bourke's description of |
| 55 | // the CONREC algorithm (a triangle with three vertices the lie |
| 56 | // on the contour level). He says, "... case g above has no really |
| 57 | // satisfactory solution and fortunately will occur rarely with |
| 58 | // real arithmetic." and then goes on to show the following image: |
| 59 | // |
| 60 | // http://paulbourke.net/papers/conrec/conrec3.gif |
| 61 | // |
| 62 | // which shows case g) in the set where no edge is drawn, agreeing |
| 63 | // with our axiom above. |
| 64 | // |
| 65 | // However, in the iteration over sh at conrec.c +44, a triangle |
| 66 | // with all vertices on the plane is given sh = {0,0,0,0,0} and |
| 67 | // then when the switch at conrec.c +93 happens, castab resolves |
| 68 | // that to case 3 for all values of m. |
| 69 | // |
| 70 | // This is fixed by replacing castab/cases[1][1][1] with 0. |
| 71 | cases = [3][3][3]int{ |
| 72 | {{0, 0, 8}, {0, 2, 5}, {7, 6, 9}}, |
| 73 | {{0, 3, 4}, {1, 0, 1}, {4, 3, 0}}, |
| 74 | {{9, 6, 7}, {5, 2, 0}, {8, 0, 0}}, |
| 75 | } |
| 76 | ) |
| 77 | |
| 78 | c, r := g.Dims() |
| 79 | for i := range c - 1 { |
| 80 | for j := range r - 1 { |
| 81 | dmin := math.Min( |
| 82 | math.Min(g.Z(i, j), g.Z(i, j+1)), |
| 83 | math.Min(g.Z(i+1, j), g.Z(i+1, j+1)), |
| 84 | ) |
| 85 | |
| 86 | dmax := math.Max( |
| 87 | math.Max(g.Z(i, j), g.Z(i, j+1)), |
| 88 | math.Max(g.Z(i+1, j), g.Z(i+1, j+1)), |
| 89 | ) |
| 90 | |
| 91 | if dmax < heights[0] || heights[len(heights)-1] < dmin { |
| 92 | continue |