class GSL::Contour

A class for making contours of a function on a regular two dimensional grid. If contours of scattered data are required, see GSL::ScatterInterp#to_contour.

Constants

VALID_CONNECTIONS

Edges: 1 __

0 | 4 \ /5  | 2
  | 7 / \ 6 |
    3 __

Attributes

keep_path_data[RW]

Public Class Methods

alloc(x, y, grid) click to toggle source

Create a new Contour object. x and y are vectors of coordinates, and grid is a matrix of values on those coordinates.

# File lib/gsl_extras.rb, line 279
def self.alloc(x, y, grid)
        new(x, y, grid)
end
new(x, y, grid) click to toggle source
# File lib/gsl_extras.rb, line 285
        def initialize(x, y, grid)
                @x = x; @y=y; @grid=grid
#                       p @grid, @x, @y
                raise ArgumentError.new("Unmatching data sizes: #{x.size}, #{y.size}, #{grid.shape}") unless [x.size, y.size] == grid.shape
                @adaptive = false                     
        end

Public Instance Methods

contours(*values) click to toggle source

Create a series of contours at the given values. Returns a hash of {value => array_of_contours}. The array_of_contours is a an array of arrays, where each array is a list of [x, y] coordinates along the contour.

# File lib/gsl_extras.rb, line 309
        def contours(*values)
                (values = (0..values[0]+1).to_a.map{|i| i.to_f * (@grid.max - @grid.min) / ( values[0]+1) + @grid.min}; values.pop; values.shift) if values.size==1 and values[0].kind_of? Integer
                cons = values.inject({}){|hash, val| hash[val] = []; hash}
#                       p cons
                for i in 0...((@x.size / 2.0).ceil - 1)
                        for j in 0...((@y.size / 2.0).ceil - 1)
                                analyse_cell(i*2, j*2, cons)
                        end
                end
#                       pp cons
                cons.keys.each{|val| cons[val] = connect_contours(cons[val])}
                @last_contours = cons
        end
graphkit(*args) click to toggle source

Create a GraphKit object of the contours.

# File lib/gsl_extras.rb, line 325
def graphkit(*args)
        if args.size == 0
                conts = @last_contours
        else
                conts = contours(*args)
        end
        graphs = conts.map do |val, cons|
                unless cons[0]
                        nil
                else
                        (cons.map do |con|
#                              p con
                                contour = con.transpose
                                kit = CodeRunner::GraphKit.autocreate({x: {data: contour[0]}, y: {data: contour[1], title: val.to_s}})
                                kit.data[0].with = "l"
                                kit
                        end).sum
                end
        end
        graphs.compact.reverse.sum
end

Private Instance Methods

analyse_cell(i, j, cons) click to toggle source
# File lib/gsl_extras.rb, line 387
        def analyse_cell(i, j, cons)
                crossed_edges = get_crossed_edges(i, j, cons)
                if @adaptive and crossed_edges.values.find_all{|crossings| crossings.size>0}.size > 0
                        @adaptive_matrix ||= GSL::Matrix.alloc(@adaption_scale+1, @adaption_scale+1)
                        @adaptive_x ||= GSL::Vector.alloc(@adaption_scale+1)
                        @adaptive_y ||= GSL::Vector.alloc(@adaption_scale+1)
                        dx = (@x[i+2] - @x[i])/@adaption_scale.to_f
                        dy = (@y[j+2] - @y[j])/@adaption_scale.to_f
                        x = @x[i]
                        y = @y[j]
                        for ii in 0...@adaption_scale+1
                                for jj in 0...@adaption_scale+1
                                        @adaptive_x[ii] = x + ii * dx
                                        @adaptive_y[jj] = y + jj * dy
                                        @adaptive_matrix[ii, jj] = @func.eval(@adaptive_x[ii], @adaptive_y[jj])
                                end
                        end
                        cell_contour = Contour.alloc(@adaptive_x, @adaptive_y, @adaptive_matrix)
                        cell_contour.keep_path_data = true
                        if @multi_adaptive and @next_adaption_scale > 1
                                #p @next_adaption_scale
                                #p "#{@adaptive_x.max}, #{@adaptive_x.min}, #{@x[i]}, #{@x[i+2]}"
                                cell_contour.set_adaptive(@func, @next_adaption_scale, true)
                        end

                        cell_cons = cell_contour.contours(*cons.keys)
                        #kit = cell_contour.graphkit
                        #kit.gnuplot(xrange: [@x[i], @x[i+2]], yrange: [@y[j], @y[j+2]], style: "data lp")
                        #p cell_cons
                        cell_cons.each do |val, cell_val_contours|
                                cell_val_contours.each do |path|
                                        # we have to relabel the path from the fine scale contour so that it fits with the course scale labelling system. only the beginning and end of the path matters
                                        path_start = path[0][0]  
                                        if path_start[0] == @x[i]
                                                        path[0][1] = [i, j, 1]
                                        elsif path_start[0] == @x[i+2]
                                                  path[0][1] = [i, j, 3]
                                        else
                                                if path_start[1] == @y[j]
                                                  path[0][1] = [i, j, 0]
                                                elsif path_start[1] == @y[j+2]
                                                  path[0][1] = [i, j, 2]
                                                else
                                                        raise "Could not find path_start; #{path_start}; x: #{@x[i]}, #{@x[i+2]}; y: #{@y[j]}, #{@y[j+2]} "
                                                end
                                        end
                                        path_end = path[-1][0]  
                                        if path_end[0] == @x[i]
                                                        path[-1][1] = [i, j, 1]
                                        elsif path_end[0] == @x[i+2]
                                                  path[-1][1] = [i, j, 3]
                                        else
                                                if path_end[1] == @y[j]
                                                  path[-1][1] = [i, j, 0]
                                                elsif path_end[1] == @y[j+2]
                                                  path[-1][1] = [i, j, 2]
                                                else
                                                        raise "Could not find path_end #{path_end}; x: #{@x[i]}, #{@x[i+2]}; y: #{@y[j]}, #{@y[j+2]} "
                                                end
                                        end

                                        cons[val].push path
                                end
                        end
                        #kit.close
                else
                        crossed_edges.each do |val, crossings|
                                outer = crossings.keys.find_all{|edge| edge < 4}
                                inner = crossings.keys.find_all{|edge| edge > 4}
                                next if outer.size == 0 and inner.size == 0
                                VALID_CONNECTIONS.each do |connection|
                                        path = crossings.values_at(*connection).compact.zip(connection.map{|edge| [i, j, edge]})
        #                                      p path
                                        next if path.size != connection.size
                                        connection.each{|edge| crossings.delete(edge)}
                                        cons[val].push path
                                        
                                end
                        end

#                               p val, crossings, cons
                end
        end
connect_contours(contours) click to toggle source
# File lib/gsl_extras.rb, line 473
        def connect_contours(contours)
                return contours if contours.size == 1
                loop do
                        catch(:restart) do 
#                               old_contours = contours
#                               contours = []
                        joined = []
                        for i in 0...contours.size
                                break if i >= contours.size
                                for j in i+1...contours.size
                                        break if j >= contours.size
                                        coni = contours[i]
                                        conj = contours[j]
                                        if joins?(coni[-1], conj[0])
                                                contours[i] = coni + conj
                                                contours.delete_at(j)
#                                                       redo
                                                throw(:restart)
                                        elsif joins?(coni[0], conj[-1])
                                                contours[i] = conj + coni
                                                contours.delete_at(j)
#                                                       redo
                                                throw(:restart)
                                        elsif joins?(coni[0], conj[0])
                                                contours[i] = coni.reverse + conj
                                                contours.delete_at(j)
#                                                       redo
                                                throw(:restart)
                                        elsif joins?(coni[-1], conj[-1])
                                                contours[i] = conj + coni.reverse
                                                contours.delete_at(j)
#                                                       redo
                                                throw(:restart)
                                        end


                                        
                                end
                        end
                        contours.each{|con| con.uniq!}
                        contours.map!{|con| con.map{|point| point[0]}} unless @keep_path_data
                        return contours
                        end
                        
                end
        end
get_crossed_edges(i, j, cons) click to toggle source
# File lib/gsl_extras.rb, line 355
        def get_crossed_edges(i, j, cons)
                ce = {}
                edges = {0=>[[i, j], [i+2, j]],
                                1=>[[i, j], [i, j+2]],
                                2=>[[i, j+2], [i+2, j+2]],
                                3=>[[i+2, j], [i+2, j+2]],
                                4=>[[i, j], [i+1, j+1]],
                                5=>[[i, j+2], [i+1, j+1]],
                                6=>[[i+2, j+2], [i+1, j+1]],
                                7=>[[i+2, j], [i+1, j+1]]
                        }    
                cons.keys.each do |val|
                        ce[val] = {}
                        edges.each do |edge, (start, fin)|
                                bounds = [@grid[*start], @grid[*fin]]
#                                       p edge, bounds if bounds.max < 4
                                if val <= bounds.max and val > bounds.min
                                        dx = @x[fin[0]] - @x[start[0]]
                                        dy = @y[fin[1]] - @y[start[1]]
                                        df = bounds[1] - bounds[0]
                                        xcross = @x[start[0]] + (val-bounds[0]) / df * dx
                                        ycross = @y[start[1]] + (val-bounds[0]) / df * dy
                                        ce[val][edge] = [xcross, ycross]
                                end
                                
                        end
                end
                ce
        end
joins?(path, contour) click to toggle source
# File lib/gsl_extras.rb, line 524
        def joins?(path, contour)      
#                       p @x.size
#                       dx = (@x.subvector(1, @x.size - 1) - @x.subvector(0, @x.size - 1)).sum / @x.size #The average dx between gridpoints
#                       dy = (@y.subvector(1, @y.size - 1) - @y.subvector(0, @y.size - 1)).sum / @y.size
#                       eps = 1.0
# #                     return true
#                       return ((path[0] - contour[0]).abs < dx.abs * eps and (path[1] - contour[1]).abs < dy.abs * eps)
#                       p dx; exit
                pi, pj, pedge = path[1]
                ci, cj, cedge = contour[1]
                joins = false
                joins = ((pi == ci and pj == cj + 2 and pedge == 0 and cedge == 2) or
                        (pi == ci and pj == cj - 2 and pedge == 2 and cedge == 0) or
                        (pi == ci + 2 and pj == cj  and pedge == 1 and cedge == 3) or
                        (pi == ci - 2 and pj == cj  and pedge == 3 and cedge == 1))
#                       unless joins
#                               p path[1], contour[1]
#                               STDIN.gets
#                       end
#                       if path[1] == [2,10,1]
#                               p path[1], contour[1], joins
#                               STDIN.gets
#                       end
                return joins
        end