Switch to ndarray and collapse layers

src/app.coffee

@@ -1,4 +1,5 @@
-
+ndarray = require('ndarray')
+ops = require('ndarray-ops')
 window.Viewer or= {}
 
 ### VARIOUS HELPFUL FUNCTIONS ###
@@ -56,12 +57,12 @@ window.Viewer = class Viewer
       v.clear()
       # Paint all layers. Note the reversal of layer order to ensure 
       # top layers get painted last.
-      for l in @layerList.layers.slice(0).reverse()
-        v.paint(l) if l.visible
+      v.paint(@layerList.layers.slice(0).reverse())
+
       # v.paint((l for l in @layerList.layers.slice(0).reverse() if l.visible))
       v.drawCrosshairs()
       v.drawLabels()
-    $(@).trigger("beforePaint")
+    $(@).trigger("afterPaint")
     return true
 
 
@@ -316,9 +317,9 @@ window.Viewer = class Viewer
       [x, y, z] = coords
     else
       [x, y, z] = @coords_ijk
-    return (l.image.data[z][y][x] for l in @layerList.layers) if all
+    return (l.image.data.get(z, y, x) for l in @layerList.layers) if all
     layer = if layer? then @layerList.layers[layer] else @layerList.activeLayer
-    layer.image.data[z][y][x]
+    layer.image.data.get(z, y, x)
 
   # Takes dimension and x/y as input and returns x/y/z viewer coordinates
   viewer2dTo3d: (dim, cx, cy = null) ->

src/models.coffee

@@ -9,31 +9,21 @@ class Image
     # Images loaded from a binary volume already have 3D data, and we 
     # just need to clean up values and swap axes (to reverse x and z 
     # relative to xtk).
+    @data = ndarray(new Float32Array(@x*@y*@z), [@x, @y, @z])
+
     if 'data3d' of data
-      @min = 0
-      @max = 0
-      @data = []
       for i in [0...@x]
-        @data[i] = []
         for j in [0...@y]
-          @data[i][j] = []
           for k in [0...@z]
-            value = Math.round(data.data3d[i][j][k]*100)/100
-            @max = value if value > @max
-            @min = value if value < @min
-            @data[i][j][k] = value
+            @data.set(i, j, k, data.data3d[i][j][k])
 
     # Load from JSON format. The format is kind of clunky and could be improved.
     else if 'values' of data
-      [@max, @min] = [data.max, data.min]
       vec = Transform.jsonToVector(data)
-      @data = Transform.vectorToVolume(vec, [@x, @y, @z])
+      @data = ndarray(vec, [@x, @y, @z])
 
-    # Otherwise initialize a blank image.
-    else
-      @min = 0
-      @max = 0
-      @data = @empty()
+    @min = ops.inf(@data)
+    @max = ops.sup(@data)
 
     # If peaks are passed, construct spheres around them
     if 'peaks' of data
@@ -42,18 +32,6 @@ class Image
              # setting to twice the value seems to work.
 
 
-  # Return an empty volume of current image dimensions
-  empty: () ->
-    vol = []
-    for i in [0...@x]
-      vol[i] = []
-      for j in [0...@y]
-        vol[i][j] = []
-        for k in [0...@z]
-          vol[i][j][k] = 0
-    return vol
-
-
   # Add a sphere of radius r at the provided coordinates. Coordinates are specified
   # in image space (i.e., where x/y/z are indexed from 0 to the number of voxels in
   # each plane).
@@ -68,37 +46,29 @@ class Image
         for k in [-r..r]
           continue if (z-k) < 0 or (z+k) > (@z - 1)
           dist = i*i + j*j + k*k
-          @data[i+x][j+y][k+z] = value if dist < r*r
+          @data.set(i+x, j+y, k+z, value) if dist < r*r
     return false
 
 
   # Need to implement resampling to allow display of images of different resolutions
   resample: (newx, newy, newz) ->
 
-
   # Slice the volume along the specified dimension (0 = x, 1 = y, 2 = z) at the
   # specified index and return a 2D array.
   slice: (dim, index) ->
     switch dim
-        when 0
-          slice = []
-          for i in [0...@x]
-            slice[i] = []
-            for j in [0...@y]
-              slice[i][j] = @data[i][j][index]
-        when 1
-          slice = []
-          for i in [0...@x]
-            slice[i] = @data[i][index]
-        when 2
-          slice = @data[index]
+      when 0
+        slice = @data.pick(null, null, index)
+      when 1
+        slice = @data.pick(null, index, null)
+      when 2
+        slice = @data.pick(index, null, null)
     return slice
 
   dims: ->
     return [@x, @y, @z]
 
 
-
 class Layer
 
   # In addition to basic properties we attach to current Layer instance,
@@ -123,34 +93,25 @@ class Layer
 
     @name = options.name
     @sign = options.sign
-    @colorMap = @setColorMap(options.colorPalette)
     @visible = options.visible
     @threshold = @setThreshold(options.negativeThreshold, options.positiveThreshold)
     @opacity = options.opacity
     @download = options.download
     @intent = options.intent
     @description = options.description
-
+    @setColorMap(options.colorPalette)
 
   hide: ->
     @visible = false
 
-
   show: ->
     @visible = true
 
-
   toggle: ->
     @visible = !@visible
 
-
-  slice: (view, viewer) ->
-    # get the right 2D slice from the Image
-    data = @image.slice(view.dim, viewer.coords_ijk[view.dim])
-    # Threshold if needed
-    data = @threshold.mask(data)
-    return data
-
+  slice: (dim) ->
+    @image.slice(dim, viewer.coords_ijk[dim])
 
   setColorMap: (palette = null, steps = null) ->
     @palette = palette
@@ -174,28 +135,29 @@ class Layer
       min = if @sign == 'positive' then 0 else @image.min
       max = if @sign == 'negative' then 0 else @image.max
     @colorMap = new ColorMap(min, max, palette, steps)
+    # @colorData = @colorMap.mapVolume(@image.data)
 
 
   setThreshold: (negThresh = 0, posThresh = 0) ->
     @threshold = new Threshold(negThresh, posThresh, @sign)
 
-
   # Update the layer's settings from provided object.
   update: (settings) ->
+
     # Handle settings that take precedence first
     @sign = settings['sign'] if 'sign' of settings
 
-    # Now everything else
+    # Now everything else--ignoring settings that haven't changed
     nt = 0
     pt = 0
     for k, v of settings
       switch k
-        when 'colorPalette' then @setColorMap(v)
-        when 'opacity' then @opacity = v
-        when 'image-intent' then @intent = v
+        when 'colorPalette' then @setColorMap(v) if @palette != v
+        when 'opacity' then @opacity = v if @opacity != v
+        when 'image-intent' then @intent = v if @intent != v
         when 'pos-threshold' then pt = v
         when 'neg-threshold' then nt = v
-        when 'description' then @description = v
+        when 'description' then @description = v if @description != v
     @setThreshold(nt, pt, @sign)
 
 
@@ -311,6 +273,70 @@ class LayerList
     @layers = newLayers
 
 
+class ColorMap
+
+  @hexToRgb = (hex) ->
+    result = /^#?([a-f\d]{2})([a-f\d]{2})([a-f\d]{2})$/i.exec(hex)
+    if result
+      [parseInt(result[1], 16),
+      parseInt(result[2], 16),
+      parseInt(result[3], 16)]
+    else [NaN, NaN, NaN]
+
+  @componentToHex = (c) ->
+    hex = c.toString(16)
+    (if hex.length is 1 then "0" + hex else hex)
+
+  @rgbToHex = (rgb) ->
+    "#" + componentToHex(rgb[0]) + componentToHex(rgb[1]) + componentToHex(rgb[2])
+
+  # For now, palettes are hard-coded. Should eventually add facility for
+  # reading in additional palettes from file and/or creating them in-browser.
+  @PALETTES =
+    grayscale: ['#000000','#303030','gray','silver','white']
+  # Add monochrome palettes
+  basic = ['red', 'green', 'blue', 'yellow', 'purple', 'lime', 'aqua', 'navy']
+  for col in basic
+    @PALETTES[col] = ['black', col, 'white']
+  # Add some other palettes
+  $.extend(@PALETTES, {
+    'intense red-blue': ['#053061', '#2166AC', '#4393C3', '#F7F7F7', '#D6604D', '#B2182B', '#67001F']
+    'red-yellow-blue': ['#313695', '#4575B4', '#74ADD1', '#FFFFBF', '#F46D43', '#D73027', '#A50026']
+    'brown-teal': ['#003C30', '#01665E', '#35978F', '#F5F5F5', '#BF812D', '#8C510A', '#543005']
+  })
+
+  constructor: (@min, @max, @palette = 'hot and cold', @steps = 40) ->
+    @range = @max - @min
+    @colors = @setColors(ColorMap.PALETTES[@palette])
+
+  # Map values to colors. Currently uses a linear mapping;  could add option
+  # to use other methods.
+  map: (data) ->
+    dims = data.shape
+    dims.push(3)
+    res = ndarray(new Float32Array(dims[0] * dims[1] * 3), dims)
+    for i in [0...data.shape[0]]
+      for j in [0...data.shape[1]]
+        v = data.get(i, j)
+        if v == 0
+          rgb = [NaN, NaN, NaN]
+        else
+          val = @colors[Math.floor(((v-@min)/@range) * @steps)]
+          rgb = ColorMap.hexToRgb(val)
+        for c in [0...3]
+          res.set(i, j, c, rgb[c])
+    return res
+
+  # Takes a set of discrete color names/descriptions and remaps them to
+  # a space with @steps different colors.
+  setColors: (colors) ->
+    rainbow = new Rainbow()
+    rainbow.setNumberRange(1, @steps)
+    rainbow.setSpectrum.apply(null, colors)
+    colors = []
+    colors.push rainbow.colourAt(i) for i in [1...@steps]
+    return colors.map (c) -> "#" + c
+
 
 # Provides thresholding/masking functionality.
 class Threshold
@@ -322,10 +348,13 @@ class Threshold
   mask: (data) ->
     return data if @posThresh is 0 and @negThresh is 0 and @sign == 'both'
     # Zero out any values below threshold or with wrong sign
-    res = []
-    for i in [0...data.length]
-      res[i] = data[i].map (v) =>
-        if (@negThresh < v < @posThresh) or (v < 0 and @sign == 'positive') or (v > 0 and @sign == 'negative') then 0 else v
+    res = ndarray(new Float32Array(data.size), data.shape)
+    for i in [0...data.shape[0]]
+      for j in [0...data.shape[1]]
+        v = data.get(i, j)
+        val =
+          if (@negThresh < v < @posThresh) or (v < 0 and @sign == 'positive') or (v > 0 and @sign == 'negative') then 0 else v
+        res.set(i, j, val)
     return res
 
 
@@ -337,32 +366,14 @@ Transform =
   # Takes compressed JSON-encoded image data as input and reconstructs
   # into a dense 1D vector, indexed from 0 to the total number of voxels.
   jsonToVector: (data) ->
-    v = new Array(data.dims[0] * data.dims[1] * data.dims[2])
+    v = new Float32Array(data.dims[0] * data.dims[1] * data.dims[2])
     v[i] = 0 for i in [0...v.length]
     for i in [0...data.values.length]
       curr_inds = data.indices[i]
       for j in [0...curr_inds.length]
           v[curr_inds[j] - 1] = data.values[i]
     return(v)
 
-  # Reshape a 1D vector of all voxels into a 3D volume with specified dims.
-  vectorToVolume: (vec, dims) ->
-    vol = []
-    for i in [0...dims[0]]
-      vol[i] = []
-      for j in [0...dims[1]]
-        vol[i][j] = []
-        for k in [0...dims[2]]
-          vol[i][j][k] = 0
-          sliceSize = dims[1] * dims[2]
-    for i in [0...vec.length]
-      continue if typeof vec[i] is `undefined`
-      x = Math.floor(i / sliceSize)
-      y = Math.floor((i - (x * sliceSize)) / dims[2])
-      z = i - (x * sliceSize) - (y * dims[2])
-      vol[x][y][z] = vec[i]
-    return(vol)
-
   # Generic coordinate transformation function that takes an input
   # set of coordinates and a matrix to use in the transformation.
   # Depends on the Sylvester library.