conrec.go

Documentation: gonum.org/v1/plot/plotter

     1  // Copyright ©2015 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package plotter
     6  
     7  import (
     8  	"math"
     9  )
    10  
    11  type point struct {
    12  	X, Y float64
    13  }
    14  
    15  type line struct {
    16  	p1, p2 point
    17  }
    18  
    19  func sect(h, p [5]float64, v1, v2 int) float64 {
    20  	return (h[v2]*p[v1] - h[v1]*p[v2]) / (h[v2] - h[v1])
    21  }
    22  
    23  // conrecLine performs an operation with a line at a given height derived
    24  // from data over the 2D box interval (i, j) to (i+1, j+1).
    25  type conrecLine func(i, j int, l line, height float64)
    26  
    27  // conrec is a Go translation of the C version of CONREC by Paul Bourke:
    28  // http://paulbourke.net/papers/conrec/conrec.c
    29  //
    30  // conrec takes g, an m×n grid function, a sorted slice of contour heights
    31  // and a conrecLine function.
    32  //
    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 := 0; i < c-1; i++ {
    80  		for j := 0; j < r-1; j++ {
    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
    93  			}
    94  
    95  			for k := 0; k < len(heights); k++ {
    96  				if heights[k] < dmin || dmax < heights[k] {
    97  					continue
    98  				}
    99  				for m := 4; m >= 0; m-- {
   100  					if m > 0 {
   101  						h[m] = g.Z(i+im[m-1], j+jm[m-1]) - heights[k]
   102  						xh[m] = g.X(i + im[m-1])
   103  						yh[m] = g.Y(j + jm[m-1])
   104  					} else {
   105  						h[0] = 0.25 * (h[1] + h[2] + h[3] + h[4])
   106  						xh[0] = 0.50 * (g.X(i) + g.X(i+1))
   107  						yh[0] = 0.50 * (g.Y(j) + g.Y(j+1))
   108  					}
   109  					switch {
   110  					case h[m] > 0:
   111  						sh[m] = 1
   112  					case h[m] < 0:
   113  						sh[m] = -1
   114  					default:
   115  						sh[m] = 0
   116  					}
   117  				}
   118  
   119  				/*
   120  				   Note: at this stage the relative heights of the corners and the
   121  				   centre are in the h array, and the corresponding coordinates are
   122  				   in the xh and yh arrays. The centre of the box is indexed by 0
   123  				   and the 4 corners by 1 to 4 as shown below.
   124  				   Each triangle is then indexed by the parameter m, and the 3
   125  				   vertices of each triangle are indexed by parameters m1,m2,and m3.
   126  				   It is assumed that the centre of the box is always vertex 2
   127  				   though this isimportant only when all 3 vertices lie exactly on
   128  				   the same contour level, in which case only the side of the box
   129  				   is drawn.
   130  
   131  				      vertex 4 +-------------------+ vertex 3
   132  				               | \               / |
   133  				               |   \    m-3    /   |
   134  				               |     \       /     |
   135  				               |       \   /       |
   136  				               |  m=2    X   m=2   |       the centre is vertex 0
   137  				               |       /   \       |
   138  				               |     /       \     |
   139  				               |   /    m=1    \   |
   140  				               | /               \ |
   141  				      vertex 1 +-------------------+ vertex 2
   142  				*/
   143  
   144  				// Scan each triangle in the box.
   145  				for m := 1; m <= 4; m++ {
   146  					m1 := m
   147  					const m2 = 0
   148  					var m3 int
   149  					if m != 4 {
   150  						m3 = m + 1
   151  					} else {
   152  						m3 = 1
   153  					}
   154  					switch cases[sh[m1]+1][sh[m2]+1][sh[m3]+1] {
   155  					case 0:
   156  						continue
   157  
   158  					case 1: // Line between vertices 1 and 2
   159  						p1 = point{X: xh[m1], Y: yh[m1]}
   160  						p2 = point{X: xh[m2], Y: yh[m2]}
   161  
   162  					case 2: // Line between vertices 2 and 3
   163  						p1 = point{X: xh[m2], Y: yh[m2]}
   164  						p2 = point{X: xh[m3], Y: yh[m3]}
   165  
   166  					case 3: // Line between vertices 3 and 1
   167  						p1 = point{X: xh[m3], Y: yh[m3]}
   168  						p2 = point{X: xh[m1], Y: yh[m1]}
   169  
   170  					case 4: // Line between vertex 1 and side 2-3
   171  						p1 = point{X: xh[m1], Y: yh[m1]}
   172  						p2 = point{X: sect(h, xh, m2, m3), Y: sect(h, yh, m2, m3)}
   173  
   174  					case 5: // Line between vertex 2 and side 3-1
   175  						p1 = point{X: xh[m2], Y: yh[m2]}
   176  						p2 = point{X: sect(h, xh, m3, m1), Y: sect(h, yh, m3, m1)}
   177  
   178  					case 6: // Line between vertex 3 and side 1-2
   179  						p1 = point{X: xh[m3], Y: yh[m3]}
   180  						p2 = point{X: sect(h, xh, m1, m2), Y: sect(h, yh, m1, m2)}
   181  
   182  					case 7: // Line between sides 1-2 and 2-3
   183  						p1 = point{X: sect(h, xh, m1, m2), Y: sect(h, yh, m1, m2)}
   184  						p2 = point{X: sect(h, xh, m2, m3), Y: sect(h, yh, m2, m3)}
   185  
   186  					case 8: // Line between sides 2-3 and 3-1
   187  						p1 = point{X: sect(h, xh, m2, m3), Y: sect(h, yh, m2, m3)}
   188  						p2 = point{X: sect(h, xh, m3, m1), Y: sect(h, yh, m3, m1)}
   189  
   190  					case 9: // Line between sides 3-1 and 1-2
   191  						p1 = point{X: sect(h, xh, m3, m1), Y: sect(h, yh, m3, m1)}
   192  						p2 = point{X: sect(h, xh, m1, m2), Y: sect(h, yh, m1, m2)}
   193  
   194  					default:
   195  						panic("cannot reach")
   196  					}
   197  
   198  					fn(i, j, line{p1: p1, p2: p2}, heights[k])
   199  				}
   200  			}
   201  		}
   202  	}
   203  }
   204
View as plain text