contour.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  	"image/color"
     9  	"math"
    10  	"sort"
    11  
    12  	"gonum.org/v1/plot"
    13  	"gonum.org/v1/plot/palette"
    14  	"gonum.org/v1/plot/vg"
    15  	"gonum.org/v1/plot/vg/draw"
    16  )
    17  
    18  // Contour implements the Plotter interface, drawing
    19  // a contour plot of the values in the GridXYZ field.
    20  type Contour struct {
    21  	GridXYZ GridXYZ
    22  
    23  	// Levels describes the contour heights to plot.
    24  	Levels []float64
    25  
    26  	// LineStyles is the set of styles for contour
    27  	// lines. Line styles are are applied to each level
    28  	// in order, modulo the length of LineStyles.
    29  	LineStyles []draw.LineStyle
    30  
    31  	// Palette is the color palette used to render
    32  	// the heat map. If Palette is nil or has no
    33  	// defined color, the Contour LineStyle color
    34  	// is used.
    35  	Palette palette.Palette
    36  
    37  	// Underflow and Overflow are colors used to draw
    38  	// contours outside the dynamic range defined
    39  	// by Min and Max.
    40  	Underflow color.Color
    41  	Overflow  color.Color
    42  
    43  	// Min and Max define the dynamic range of the
    44  	// heat map.
    45  	Min, Max float64
    46  }
    47  
    48  // NewContour creates as new contour plotter for the given data, using
    49  // the provided palette. If levels is nil, contours are generated for
    50  // the 0.01, 0.05, 0.25, 0.5, 0.75, 0.95 and 0.99 quantiles.
    51  // If g has Min and Max methods that return a float, those returned
    52  // values are used to set the respective Contour fields.
    53  // If the returned Contour is used when Min is greater than Max, the
    54  // Plot method will panic.
    55  func NewContour(g GridXYZ, levels []float64, p palette.Palette) *Contour {
    56  	var min, max float64
    57  	type minMaxer interface {
    58  		Min() float64
    59  		Max() float64
    60  	}
    61  	switch g := g.(type) {
    62  	case minMaxer:
    63  		min, max = g.Min(), g.Max()
    64  	default:
    65  		min, max = math.Inf(1), math.Inf(-1)
    66  		c, r := g.Dims()
    67  		for i := 0; i < c; i++ {
    68  			for j := 0; j < r; j++ {
    69  				v := g.Z(i, j)
    70  				if math.IsNaN(v) {
    71  					continue
    72  				}
    73  				min = math.Min(min, v)
    74  				max = math.Max(max, v)
    75  			}
    76  		}
    77  	}
    78  
    79  	if len(levels) == 0 {
    80  		levels = quantilesR7(g, defaultQuantiles)
    81  	}
    82  
    83  	return &Contour{
    84  		GridXYZ:    g,
    85  		Levels:     levels,
    86  		LineStyles: []draw.LineStyle{DefaultLineStyle},
    87  		Palette:    p,
    88  		Min:        min,
    89  		Max:        max,
    90  	}
    91  }
    92  
    93  // Default quantiles for case where levels is not explicitly set.
    94  var defaultQuantiles = []float64{0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99}
    95  
    96  // quantilesR7 returns the pth quantiles of the data in g according to the R-7 method.
    97  // http://en.wikipedia.org/wiki/Quantile#Estimating_the_quantiles_of_a_population
    98  func quantilesR7(g GridXYZ, p []float64) []float64 {
    99  	c, r := g.Dims()
   100  	data := make([]float64, 0, c*r)
   101  	for i := 0; i < c; i++ {
   102  		for j := 0; j < r; j++ {
   103  			if v := g.Z(i, j); !math.IsNaN(v) {
   104  				data = append(data, v)
   105  			}
   106  		}
   107  	}
   108  	sort.Float64s(data)
   109  	v := make([]float64, len(p))
   110  	for j, q := range p {
   111  		if q == 1 {
   112  			v[j] = data[len(data)-1]
   113  		}
   114  		h := float64(len(data)-1) * q
   115  		i := int(h)
   116  		v[j] = data[i] + (h-math.Floor(h))*(data[i+1]-data[i])
   117  	}
   118  	return v
   119  }
   120  
   121  // naive is a debugging constant. If true, Plot performs no contour path
   122  // reconstruction, instead rendering each path segment individually.
   123  const naive = false
   124  
   125  // Plot implements the Plot method of the plot.Plotter interface.
   126  func (h *Contour) Plot(c draw.Canvas, plt *plot.Plot) {
   127  	if h.Min > h.Max {
   128  		panic("contour: invalid Z range: min greater than max")
   129  	}
   130  
   131  	if naive {
   132  		h.naivePlot(c, plt)
   133  		return
   134  	}
   135  
   136  	var pal []color.Color
   137  	if h.Palette != nil {
   138  		pal = h.Palette.Colors()
   139  	}
   140  
   141  	trX, trY := plt.Transforms(&c)
   142  
   143  	// Collate contour paths and draw them.
   144  	//
   145  	// The alternative naive approach is to draw each line segment as
   146  	// conrec returns it. The integrated path approach allows graphical
   147  	// optimisations and is necessary for contour fill shading.
   148  	cp := contourPaths(h.GridXYZ, h.Levels, trX, trY)
   149  
   150  	// ps is a palette scaling factor to scale the palette uniformly
   151  	// across the given levels. This enables a discordance between the
   152  	// number of colours and the number of levels. Sorting is not
   153  	// necessary since contourPaths sorts the levels as a side effect.
   154  	ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0])
   155  	if len(h.Levels) == 1 {
   156  		ps = 0
   157  	}
   158  
   159  	for i, z := range h.Levels {
   160  		if math.IsNaN(z) {
   161  			continue
   162  		}
   163  		for _, pa := range cp[z] {
   164  			if isLoop(pa) {
   165  				pa.Close()
   166  			}
   167  
   168  			style := h.LineStyles[i%len(h.LineStyles)]
   169  			var col color.Color
   170  			switch {
   171  			case z < h.Min:
   172  				col = h.Underflow
   173  			case z > h.Max:
   174  				col = h.Overflow
   175  			case len(pal) == 0:
   176  				col = style.Color
   177  			default:
   178  				col = pal[int((z-h.Levels[0])*ps+0.5)] // Apply palette scaling.
   179  			}
   180  			if col != nil && style.Width != 0 {
   181  				c.SetLineStyle(style)
   182  				c.SetColor(col)
   183  				c.Stroke(pa)
   184  			}
   185  		}
   186  	}
   187  }
   188  
   189  // naivePlot implements a naive rendering approach for contours.
   190  // It is here as a debugging mode since it simply draws line segments
   191  // generated by conrec without further computation.
   192  func (h *Contour) naivePlot(c draw.Canvas, plt *plot.Plot) {
   193  	var pal []color.Color
   194  	if h.Palette != nil {
   195  		pal = h.Palette.Colors()
   196  	}
   197  
   198  	trX, trY := plt.Transforms(&c)
   199  
   200  	// Sort levels prior to palette scaling since we can't depend on
   201  	// sorting as a side effect from calling contourPaths.
   202  	sort.Float64s(h.Levels)
   203  	// ps is a palette scaling factor to scale the palette uniformly
   204  	// across the given levels. This enables a discordance between the
   205  	// number of colours and the number of levels.
   206  	ps := float64(len(pal)-1) / (h.Levels[len(h.Levels)-1] - h.Levels[0])
   207  	if len(h.Levels) == 1 {
   208  		ps = 0
   209  	}
   210  
   211  	levelMap := make(map[float64]int)
   212  	for i, z := range h.Levels {
   213  		levelMap[z] = i
   214  	}
   215  
   216  	// Draw each line segment as conrec generates it.
   217  	var pa vg.Path
   218  	conrec(h.GridXYZ, h.Levels, func(_, _ int, l line, z float64) {
   219  		if math.IsNaN(z) {
   220  			return
   221  		}
   222  
   223  		pa = pa[:0]
   224  
   225  		x1, y1 := trX(l.p1.X), trY(l.p1.Y)
   226  		x2, y2 := trX(l.p2.X), trY(l.p2.Y)
   227  
   228  		pt1 := vg.Point{X: x1, Y: y1}
   229  		pt2 := vg.Point{X: x2, Y: y2}
   230  		if !c.Contains(pt1) || !c.Contains(pt2) {
   231  			return
   232  		}
   233  
   234  		pa.Move(pt1)
   235  		pa.Line(pt2)
   236  		pa.Close()
   237  
   238  		style := h.LineStyles[levelMap[z]%len(h.LineStyles)]
   239  		var col color.Color
   240  		switch {
   241  		case z < h.Min:
   242  			col = h.Underflow
   243  		case z > h.Max:
   244  			col = h.Overflow
   245  		case len(pal) == 0:
   246  			col = style.Color
   247  		default:
   248  			col = pal[int((z-h.Levels[0])*ps+0.5)] // Apply palette scaling.
   249  		}
   250  		if col != nil && style.Width != 0 {
   251  			c.SetLineStyle(style)
   252  			c.SetColor(col)
   253  			c.Stroke(pa)
   254  		}
   255  	})
   256  }
   257  
   258  // DataRange implements the DataRange method
   259  // of the plot.DataRanger interface.
   260  func (h *Contour) DataRange() (xmin, xmax, ymin, ymax float64) {
   261  	c, r := h.GridXYZ.Dims()
   262  	return h.GridXYZ.X(0), h.GridXYZ.X(c - 1), h.GridXYZ.Y(0), h.GridXYZ.Y(r - 1)
   263  }
   264  
   265  // GlyphBoxes implements the GlyphBoxes method
   266  // of the plot.GlyphBoxer interface.
   267  func (h *Contour) GlyphBoxes(plt *plot.Plot) []plot.GlyphBox {
   268  	c, r := h.GridXYZ.Dims()
   269  	b := make([]plot.GlyphBox, 0, r*c)
   270  	for i := 0; i < c; i++ {
   271  		for j := 0; j < r; j++ {
   272  			b = append(b, plot.GlyphBox{
   273  				X: plt.X.Norm(h.GridXYZ.X(i)),
   274  				Y: plt.Y.Norm(h.GridXYZ.Y(j)),
   275  				Rectangle: vg.Rectangle{
   276  					Min: vg.Point{X: -2.5, Y: -2.5},
   277  					Max: vg.Point{X: +2.5, Y: +2.5},
   278  				},
   279  			})
   280  		}
   281  	}
   282  	return b
   283  }
   284  
   285  // isLoop returns true iff a vg.Path is a closed loop.
   286  func isLoop(p vg.Path) bool {
   287  	s := p[0]
   288  	e := p[len(p)-1]
   289  	return s.Pos == e.Pos
   290  }
   291  
   292  // contourPaths returns a collection of vg.Paths describing contour lines based
   293  // on the input data in m cut at the given levels. The trX and trY function
   294  // are coordinate transforms. The returned map contains slices of paths keyed
   295  // on the value of the contour level. contouPaths sorts levels ascending as a
   296  // side effect.
   297  func contourPaths(m GridXYZ, levels []float64, trX, trY func(float64) vg.Length) map[float64][]vg.Path {
   298  	sort.Float64s(levels)
   299  
   300  	ends := make(map[float64]endMap)
   301  	conts := make(contourSet)
   302  	conrec(m, levels, func(_, _ int, l line, z float64) {
   303  		paths(l, z, ends, conts)
   304  	})
   305  	ends = nil
   306  
   307  	// TODO(kortschak): Check that all non-loop paths have
   308  	// both ends at boundary. If any end is not at a boundary
   309  	// it may have a partner near by. Find this partner and join
   310  	// the two conts by merging the near by ends at the mean
   311  	// location. This operation is done level by level to ensure
   312  	// close contours of different heights are not joined.
   313  	// A partner should be a float error different end, but I
   314  	// suspect that is is possible for a bi- or higher order
   315  	// furcation so it may be that the path ends at middle node
   316  	// of another path. This needs to be investigated.
   317  
   318  	// Excise loops from crossed paths.
   319  	for c := range conts {
   320  		// Always try to do quick excision in production if possible.
   321  		c.exciseLoops(conts, true)
   322  	}
   323  
   324  	// Build vg.Paths.
   325  	paths := make(map[float64][]vg.Path)
   326  	for c := range conts {
   327  		paths[c.z] = append(paths[c.z], c.path(trX, trY))
   328  	}
   329  
   330  	return paths
   331  }
   332  
   333  // contourSet hold a working collection of contours.
   334  type contourSet map[*contour]struct{}
   335  
   336  // endMap holds a working collection of available ends.
   337  type endMap map[point]*contour
   338  
   339  // paths extends a conrecLine function to build a set of contours that represent
   340  // paths along contour lines. It is used as the engine for a closure where ends
   341  // and conts are closed around in a conrecLine function, and l and z are the
   342  // line and height values provided by conrec. At the end of a conrec call,
   343  // conts will contain a map keyed on the set of paths.
   344  func paths(l line, z float64, ends map[float64]endMap, conts contourSet) {
   345  	zEnds, ok := ends[z]
   346  	if !ok {
   347  		zEnds = make(endMap)
   348  		ends[z] = zEnds
   349  		c := newContour(l, z)
   350  		zEnds[l.p1] = c
   351  		zEnds[l.p2] = c
   352  		conts[c] = struct{}{}
   353  		return
   354  	}
   355  
   356  	c1, ok1 := zEnds[l.p1]
   357  	c2, ok2 := zEnds[l.p2]
   358  
   359  	// New segment.
   360  	if !ok1 && !ok2 {
   361  		c := newContour(l, z)
   362  		zEnds[l.p1] = c
   363  		zEnds[l.p2] = c
   364  		conts[c] = struct{}{}
   365  		return
   366  	}
   367  
   368  	if ok1 {
   369  		// Add l.p2 to end of l.p1's contour.
   370  		if !c1.extend(l, zEnds) {
   371  			panic("internal link")
   372  		}
   373  	} else if ok2 {
   374  		// Add l.p1 to end of l.p2's contour.
   375  		if !c2.extend(l, zEnds) {
   376  			panic("internal link")
   377  		}
   378  	}
   379  
   380  	if c1 == c2 {
   381  		return
   382  	}
   383  
   384  	// Join conts.
   385  	if ok1 && ok2 {
   386  		if !c1.connect(c2, zEnds) {
   387  			panic("internal link")
   388  		}
   389  		delete(conts, c2)
   390  	}
   391  }
   392  
   393  // path is a set of points forming a path.
   394  type path []point
   395  
   396  // contour holds a set of point lying sequentially along a contour line
   397  // at height z.
   398  type contour struct {
   399  	z float64
   400  
   401  	// backward and forward must each always have at least one entry.
   402  	backward path
   403  	forward  path
   404  }
   405  
   406  // newContour returns a contour starting with the end points of l for the
   407  // height z.
   408  func newContour(l line, z float64) *contour {
   409  	return &contour{z: z, forward: path{l.p1}, backward: path{l.p2}}
   410  }
   411  
   412  func (c *contour) path(trX, trY func(float64) vg.Length) vg.Path {
   413  	var pa vg.Path
   414  	p := c.front()
   415  	pa.Move(vg.Point{X: trX(p.X), Y: trY(p.Y)})
   416  	for i := len(c.backward) - 2; i >= 0; i-- {
   417  		p = c.backward[i]
   418  		pa.Line(vg.Point{X: trX(p.X), Y: trY(p.Y)})
   419  	}
   420  	for _, p := range c.forward {
   421  		pa.Line(vg.Point{X: trX(p.X), Y: trY(p.Y)})
   422  	}
   423  
   424  	return pa
   425  }
   426  
   427  // front returns the first point in the contour.
   428  func (c *contour) front() point { return c.backward[len(c.backward)-1] }
   429  
   430  // back returns the last point in the contour
   431  func (c *contour) back() point { return c.forward[len(c.forward)-1] }
   432  
   433  // extend adds the line l to the contour, updating the ends map. It returns
   434  // a boolean indicating whether the extension was successful.
   435  func (c *contour) extend(l line, ends endMap) (ok bool) {
   436  	switch c.front() {
   437  	case l.p1:
   438  		c.backward = append(c.backward, l.p2)
   439  		delete(ends, l.p1)
   440  		ends[l.p2] = c
   441  		return true
   442  	case l.p2:
   443  		c.backward = append(c.backward, l.p1)
   444  		delete(ends, l.p2)
   445  		ends[l.p1] = c
   446  		return true
   447  	}
   448  
   449  	switch c.back() {
   450  	case l.p1:
   451  		c.forward = append(c.forward, l.p2)
   452  		delete(ends, l.p1)
   453  		ends[l.p2] = c
   454  		return true
   455  	case l.p2:
   456  		c.forward = append(c.forward, l.p1)
   457  		delete(ends, l.p2)
   458  		ends[l.p1] = c
   459  		return true
   460  	}
   461  
   462  	return false
   463  }
   464  
   465  // reverse reverses the order of the point in a path and returns it.
   466  func (p path) reverse() path {
   467  	for i, j := 0, len(p)-1; i < j; i, j = i+1, j-1 {
   468  		p[i], p[j] = p[j], p[i]
   469  	}
   470  	return p
   471  }
   472  
   473  // connect connects the contour b with the receiver, updating the ends map.
   474  // It returns a boolean indicating whether the connection was successful.
   475  func (c *contour) connect(b *contour, ends endMap) (ok bool) {
   476  	switch c.front() {
   477  	case b.front():
   478  		delete(ends, c.front())
   479  		ends[b.back()] = c
   480  		c.backward = append(c.backward, b.backward.reverse()[1:]...)
   481  		c.backward = append(c.backward, b.forward...)
   482  		return true
   483  	case b.back():
   484  		delete(ends, c.front())
   485  		ends[b.front()] = c
   486  		c.backward = append(c.backward, b.forward.reverse()[1:]...)
   487  		c.backward = append(c.backward, b.backward...)
   488  		return true
   489  	}
   490  
   491  	switch c.back() {
   492  	case b.front():
   493  		delete(ends, c.back())
   494  		ends[b.back()] = c
   495  		c.forward = append(c.forward, b.backward.reverse()[1:]...)
   496  		c.forward = append(c.forward, b.forward...)
   497  		return true
   498  	case b.back():
   499  		delete(ends, c.back())
   500  		ends[b.front()] = c
   501  		c.forward = append(c.forward, b.forward.reverse()[1:]...)
   502  		c.forward = append(c.forward, b.backward...)
   503  		return true
   504  	}
   505  
   506  	return false
   507  }
   508  
   509  // exciseLoops finds loops within the contour that do not include the
   510  // start and end. Loops are removed from the contour and added to the
   511  // contour set. Loop detection is performed by Johnson's algorithm for
   512  // finding elementary cycles.
   513  func (c *contour) exciseLoops(conts contourSet, quick bool) {
   514  	if quick {
   515  		// Find cases we can guarantee don't need
   516  		// a complete analysis.
   517  		seen := make(map[point]struct{})
   518  		var crossOvers int
   519  		for _, p := range c.backward {
   520  			if _, ok := seen[p]; ok {
   521  				crossOvers++
   522  			}
   523  			seen[p] = struct{}{}
   524  		}
   525  		for _, p := range c.forward[:len(c.forward)-1] {
   526  			if _, ok := seen[p]; ok {
   527  				crossOvers++
   528  			}
   529  			seen[p] = struct{}{}
   530  		}
   531  		switch crossOvers {
   532  		case 0:
   533  			return
   534  		case 1:
   535  			c.exciseQuick(conts)
   536  			return
   537  		}
   538  	}
   539  
   540  	wp := append(c.backward.reverse(), c.forward...)
   541  	g := graphFrom(wp)
   542  	cycles := cyclesIn(g)
   543  	if len(cycles) == 0 {
   544  		// No further work to do but clean up after ourselves.
   545  		// We should not have reached here.
   546  		c.backward.reverse()
   547  		return
   548  	}
   549  	delete(conts, c)
   550  
   551  	// Put loops into the contour set.
   552  	for _, cyc := range cycles {
   553  		loop := wp.subpath(cyc)
   554  		conts[&contour{
   555  			z:        c.z,
   556  			backward: loop[:1:1],
   557  			forward:  loop[1:],
   558  		}] = struct{}{}
   559  	}
   560  
   561  	// Find non-loop paths and keep them.
   562  	g.remove(cycles)
   563  	paths := wp.linearPathsIn(g)
   564  	for _, p := range paths {
   565  		conts[&contour{
   566  			z:        c.z,
   567  			backward: p[:1:1],
   568  			forward:  p[1:],
   569  		}] = struct{}{}
   570  	}
   571  }
   572  
   573  // graphFrom returns a graph representing the point path p.
   574  func graphFrom(p path) graph {
   575  	g := make([]set, len(p))
   576  	seen := make(map[point]int)
   577  	for i, v := range p {
   578  		if _, ok := seen[v]; !ok {
   579  			seen[v] = i
   580  		}
   581  	}
   582  
   583  	for i, v := range p {
   584  		e, ok := seen[v]
   585  		if ok && g[e] == nil {
   586  			g[e] = make(set)
   587  		}
   588  		if i < len(p)-1 {
   589  			g[e][seen[p[i+1]]] = struct{}{}
   590  		}
   591  	}
   592  
   593  	return g
   594  }
   595  
   596  // subpath returns a subpath given the slice of point indices
   597  // into the path.
   598  func (p path) subpath(i []int) path {
   599  	pa := make(path, 0, len(i))
   600  	for _, n := range i {
   601  		pa = append(pa, p[n])
   602  	}
   603  	return pa
   604  }
   605  
   606  // linearPathsIn returns the linear paths in g created from p.
   607  // If g contains any cycles linearPaths will panic.
   608  func (p path) linearPathsIn(g graph) []path {
   609  	var pa []path
   610  
   611  	var u int
   612  	for u < len(g) {
   613  		for ; u < len(g) && len(g[u]) == 0; u++ {
   614  		}
   615  		if u == len(g) {
   616  			return pa
   617  		}
   618  		var curr path
   619  		for {
   620  			if len(g[u]) == 0 {
   621  				curr = append(curr, p[u])
   622  				pa = append(pa, curr)
   623  				if u == len(g)-1 {
   624  					return pa
   625  				}
   626  				break
   627  			}
   628  			if len(g[u]) > 1 {
   629  				panic("contour: not a linear path")
   630  			}
   631  			for v := range g[u] {
   632  				curr = append(curr, p[u])
   633  				u = v
   634  				break
   635  			}
   636  		}
   637  	}
   638  
   639  	return pa
   640  }
   641  
   642  // exciseQuick is a heuristic approach to loop excision. It does not
   643  // correctly identify loops in all cases, but those cases are likely
   644  // to be rare.
   645  func (c *contour) exciseQuick(conts contourSet) {
   646  	wp := append(c.backward.reverse(), c.forward...)
   647  	seen := make(map[point]int)
   648  	for j := 0; j < len(wp); {
   649  		p := wp[j]
   650  		if i, ok := seen[p]; ok && p != wp[0] && p != wp[len(wp)-1] {
   651  			conts[&contour{
   652  				z:        c.z,
   653  				backward: path{wp[i]},
   654  				forward:  append(path(nil), wp[i+1:j+1]...),
   655  			}] = struct{}{}
   656  			wp = append(wp[:i], wp[j:]...)
   657  			j = i + 1
   658  		} else {
   659  			seen[p] = j
   660  			j++
   661  		}
   662  	}
   663  	c.backward = c.backward[:1]
   664  	c.backward[0] = wp[0]
   665  	c.forward = wp[1:]
   666  }
   667
View as plain text