Module: FindBeads

Defined in:

lib/find_beads.rb,
lib/find_beads/version.rb,
lib/find_beads/bead_clumps.rb

Overview

– /* ***** BEGIN LICENSE BLOCK *****

* 
* Copyright (c) 2013 Colin J. Fuller
* 
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the Software), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
* 
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
* 
* THE SOFTWARE IS PROVIDED AS IS, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
* SOFTWARE.
* 
* ***** END LICENSE BLOCK ***** */

++

Defined Under Namespace

Modules: BeadClumping

Constant Summary

DEFAULT_SEG_CH =

2

DEFAULT_SEG_PL =

8

DEFAULT_BEAD_RADIUS =

24.0

DEFAULT_THREADS =

1

VERSION =

'0.9.5'

Class Method Summary

• .calculate_max_size_from_radius(rad) ⇒ Object

  Caclulates the maximum allowed size of a bead from the supplied radius.

• .calculate_min_size_from_radius(rad) ⇒ Object

  Calculates the minimum allowed size of a bead from the supplied radius.

• .centroids(mask) ⇒ Hash

  Finds the centroid of each unique-greylevel region in a mask.

• .is_on_voronoi_border?(points, ic) ⇒ Boolean

  Checks if a given coordinate would be approximately on the boundary between two regions of a Voronoi diagram of constructed from a set of points.

• .mask_from_image(im, opts) ⇒ Object

  Generates a segmented mask of beads from an image.

• .process_file(fn, opts = nil) ⇒ Object

  Processes a single file, which consists of creating a mask, quantifying regions, and writing output.

• .project_point_onto_vector(origin, point_to_project, point_on_line) ⇒ Array

  Projects a point in space onto a specified line.

• .recursive_threshold(p, im, mask) ⇒ void

  Recursively thresholds regions in a supplied mask and image using the method described in Xiong et al.

• .run_find_beads ⇒ Object

  Runs the bead finding on a file or directory, and grabs options from the command line.

• .write_output(fn_orig, quant_str, mask) ⇒ Object

  Writes the output data and mask to files.

Class Method Details

.calculate_max_size_from_radius(rad) ⇒ Object

Caclulates the maximum allowed size of a bead from the supplied radius. This is set to be slightly larger than a circle of that radius.

Parameters:

• rad (Fixnum) —

  the radius of the bead in units of pixels.

# File 'lib/find_beads.rb', line 278

def self.calculate_max_size_from_radius(rad)

  ((rad+1)**2 * 3.2).to_i

end

.calculate_min_size_from_radius(rad) ⇒ Object

Calculates the minimum allowed size of a bead from the supplied radius. This is set to be slightly smaller than a third of a circle of that radius. (Note that this is smaller than any of the returned regions should be, but making the cutoff this small is useful for dividing up clumps of beads where several rounds or recursive thresholding may make the regions quite small temporarily.)

Parameters:

• @see —

  #calculate_max_size_from_radius

# File 'lib/find_beads.rb', line 293

def self.calculate_min_size_from_radius(rad)

  (0.96* (rad+1)**2).to_i

end

.centroids(mask) ⇒ Hash

Finds the centroid of each unique-greylevel region in a mask.

Parameters:

• mask (Image) —

  the mask in which the regions appear. 0 denotes background and will not be counted.

Returns:

• (Hash) —

  a hash where keys are the greylevels of each region in the mask, and the values are two-element arrays containing the x,y-coordinates of the centroids of these regions.

# File 'lib/find_beads.rb', line 68

def self.centroids(mask)

  cens = {}

  mask.each do |ic|

    next unless mask[ic] > 0
    
    cens[mask[ic]] = [0.0, 0.0] unless cens[mask[ic]]

    cens[mask[ic]][0] += ic[:x]
    cens[mask[ic]][1] += ic[:y]
    
  end

  h = Histogram.new(mask)

  cens.each_key do |k|

    cens[k].map! { |e| e / h.getCounts(k) }

  end

  cens

end

.is_on_voronoi_border?(points, ic) ⇒ Boolean

Checks if a given coordinate would be approximately on the boundary between two regions of a Voronoi diagram of constructed from a set of points. The approximation is calculated such that no two regions in the Voronoi diagram would be 8-connected.

Parameters:

• points (Array) —

  an array of two-element arrays containing the x,y coordinates of the points on which the diagram is calculated.

• ic (ImageCoordinate) —

  an ImageCoordinate specifying the location to check.

Returns:

• (Boolean) —

  whether the point is on the border between two regions.

# File 'lib/find_beads.rb', line 134

def self.is_on_voronoi_border?(points, ic)

  x = ic[:x]
  y = ic[:y]

  closest_index = 0
  next_index = 0
  closest_dist = Float::MAX
  next_dist = Float::MAX
  

  points.each_with_index do |p, i|

    dist = Math.hypot(p[0] - x, p[1] - y)

    if dist < closest_dist then

      next_dist = closest_dist

      next_index = closest_index

      closest_dist = dist

      closest_index = i

    elsif dist < next_dist then

      next_dist = dist

      next_index = i

    end

  end

  proj_point = project_point_onto_vector(points[closest_index], [x,y], points[next_index])

  next_dist_proj = Math.hypot(points[next_index][0]-proj_point[0], points[next_index][1]-proj_point[1])
  closest_dist_proj = Math.hypot(points[closest_index][0]-proj_point[0], points[closest_index][1]-proj_point[1])


  cutoff = 1.01*Math.sqrt(2)

  if next_dist_proj - closest_dist_proj < cutoff then

    true

  else 

    false

  end

end

.mask_from_image(im, opts) ⇒ Object

Generates a segmented mask of beads from an image.

Parameters:

• im (Image) —

  the image to segment

• opts (Hash) —

  a hash of commandline arguments.

# File 'lib/find_beads.rb', line 305

def self.mask_from_image(im, opts)

  seg_ch = nil
  seg_pl = nil
  rad = nil

  if opts then
    seg_ch = opts[:segchannel]
    seg_pl = opts[:segplane]
    rad = opts[:beadradius]
  else
    seg_ch = DEFAULT_SEG_CH
    seg_pl = DEFAULT_SEG_PL
    rad = DEFAULT_BEAD_RADIUS
  end
  
  min_size = calculate_min_size_from_radius(rad)
  max_size = calculate_max_size_from_radius(rad)
  
  sizes = ImageCoordinate.cloneCoord(im.getDimensionSizes)

  sizes[:c] = 1
  sizes[:z] = 1

  im0 = ImageCoordinate.createCoordXYZCT(0,0,0,0,0)

  im0[:c] = seg_ch
  im0[:z] = seg_pl

  to_seg = im.subImage(sizes, im0).writableInstance

  p = RImageAnalysisTools.create_parameter_dictionary(min_size: min_size, max_size: max_size)

  im_cp = ImageFactory.create_writable(to_seg)

  mstf = MaximumSeparabilityThresholdingFilter.new

  lf = LabelFilter.new

  saf = SizeAbsoluteFilter.new

  filters = [] 

  filters << mstf

  filters << lf

  filters.each do |f|

    f.setParameters(p)
    f.setReferenceImage(im_cp)
    f.apply(to_seg)

  end

  recursive_threshold(p, im_cp, to_seg)

  saf.setParameters(p)
  
  saf.apply(to_seg)

  cens = centroids(to_seg)

  final_mask = ImageFactory.create_writable(to_seg)

  radius = rad

  final_mask.each do |ic|

    final_mask[ic] = 0
    
  end

  final_mask.each do |ic|

    x = ic[:x]
    y = ic[:y]

    cens.each_key do |k|

      if Math.hypot(cens[k][0] - x, cens[k][1] - y) <= radius then

        final_mask[ic] = k

      end

    end

  end


  final_mask.each do |ic|

    next unless final_mask[ic] > 0

    if is_on_voronoi_border?(cens.values, ic) then

      final_mask[ic] = 0

    end

  end

  lf.apply(final_mask)
  
  saf.apply(final_mask)

  lf.apply(final_mask)

  final_mask

end

.process_file(fn, opts = nil) ⇒ Object

Processes a single file, which consists of creating a mask, quantifying regions, and writing output.

Parameters:

• fn (String) —

  the filename of the image to process

• opts (Hash) (defaults to: nil) —

  a hash of command line options.

# File 'lib/find_beads.rb', line 464

def self.process_file(fn, opts=nil)

  puts "processing #{fn}"

  im = RImageAnalysisTools.get_image(fn)
  
  mask = mask_from_image(im, opts)

  proj = Java::edu.stanford.cfuller.imageanalysistools.frontend.MaximumIntensityProjection.projectImage(im)

  ims = proj.splitChannels

  is = ImageSet.new(ParameterDictionary.emptyDictionary)

  ims.each do |imc|

    is.addImageWithImage(imc)

  end

  met = IntensityPerPixelMetric.new

  q = met.quantify(mask, is)

  outdat = Java::edu.stanford.cfuller.imageanalysistools.frontend.LocalAnalysis.generateDataOutputString(q, nil)

  write_output(fn, outdat, mask)

end

.project_point_onto_vector(origin, point_to_project, point_on_line) ⇒ Array

Projects a point in space onto a specified line.

Parameters:

• origin (Array) —

  the point in space the will serve as the origin for the purposes of the projection (this should be on the line)

• point_to_project (Array) —

  the point in space that will be projected this should be in the same coordinate system in which the origin is specified, not relative to the origin

• point_on_line (Array) —

  another point on the line specified in the same coordinate system in which the origin is specified, not relative to the origin

Returns:

• (Array) —

  the projected point (in the same coordinate system as the other points were specified)

# File 'lib/find_beads.rb', line 111

def self.project_point_onto_vector(origin, point_to_project, point_on_line)

  unit_vec = (Vector[*point_on_line] - Vector[*origin]).normalize

  proj_vec = Vector[*point_to_project] - Vector[*origin]

  projected = unit_vec * (proj_vec.inner_product(unit_vec)) + Vector[*origin]

  projected.to_a

end

.recursive_threshold(p, im, mask) ⇒ void

This method returns an undefined value.

Recursively thresholds regions in a supplied mask and image using the method described in Xiong et al. (DOI: 10.1109/ICIP.2006.312365).

Parameters:

• p (ParameterDictionary) —

  a ParameterDictionary specifying max_size and min_size parameters, which control the maximum size of regions before they are recursively thresholded to break them up, and the minimum size of regions before they are discarded.

• im (Image) —

  the original image being segmented. This will not be modified.

• mask (Image) —

  the initial segmentation mask for the supplied image. Regions in this mask may be divided up or discarded.

# File 'lib/find_beads.rb', line 201

def self.recursive_threshold(p, im, mask)

  h = Histogram.new(mask)

  changed = false

  discard_list = {}

  1.upto(h.getMaxValue) do |i|

    if h.getCounts(i) > p[:max_size].to_i then

      values = []

      im.each do |ic|

        if mask[ic] == i then

          values << im[ic]

        end

      end

      thresh = RImageAnalysisTools.graythresh(values)

      im.each do |ic|

        if mask[ic] == i and im[ic] <= thresh then

          mask[ic] = 0

          changed = true

        end

      end

    elsif h.getCounts(i) > 0 and h.getCounts(i) < p[:min_size].to_i then

      discard_list[i] = true

    end

  end

  im.each do |ic|

    if discard_list[im[ic]] then

      mask[ic] = 0

      changed = true
      
    end

  end


  if changed then

    lf = LabelFilter.new

    lf.apply(mask)

    recursive_threshold(p, im, mask)

  end

end

.run_find_beads ⇒ Object

Runs the bead finding on a file or directory, and grabs options from the command line.

# File 'lib/find_beads.rb', line 497

def self.run_find_beads

  opts = Trollop::options do

    opt :dir, "Directory to process", :type => :string
    opt :file, "File to process", :type => :string
    opt :segchannel, "Channel on which to segment (0-indexed)", :type => :integer, :default => DEFAULT_SEG_CH
    opt :segplane, "Plane on which to segment (0-indexed)", :type => :integer, :default => DEFAULT_SEG_PL
    opt :max_threads, "Maximum number of paralell execution threads", :type => :integer, :default => DEFAULT_THREADS
    opt :beadradius, "Radius of the bead in pixels", :type => :float, :default => DEFAULT_BEAD_RADIUS

  end

  if opts[:dir] then

    fod = opts[:dir]

    sleep_time_s = 0.5
    threads = []

    Dir.foreach(fod) do |f|

      until threads.count { |t| t.alive? } < opts[:max_threads] do
        sleep sleep_time_s
      end

      fn = File.expand_path(f, fod)

      if File.file?(fn) then

        begin

          threads << Thread.new do 

            process_file(fn, opts)

          end

        rescue Exception => e

          puts "Unable to process #{fn}:"
          puts e.message

        end
        
      end

    end

    threads.each { |t| t.join }

  end



  if opts[:file] then

    process_file(opts[:file], opts)
    
  end

end

.write_output(fn_orig, quant_str, mask) ⇒ Object

Writes the output data and mask to files.

Parameters:

• fn_orig (String) —

  the original filename of the image being segmented/quantified.

• quant_str (String) —

  the quantification data to write

• mask (Image) —

  the mask to write

# File 'lib/find_beads.rb', line 425

def self.write_output(fn_orig, quant_str, mask)

  mask_dir = "output_mask"

  quant_dir = "quantification"

  mask_ext = "_mask.ome.tif"

  quant_ext = "_quant.txt"

  dir = File.dirname(fn_orig)

  base = File.basename(fn_orig)

  base = base.gsub(".ome.tif", "")

  mask_dir = File.expand_path(mask_dir, dir)

  quant_dir = File.expand_path(quant_dir, dir)

  Dir.mkdir(mask_dir) unless Dir.exist?(mask_dir)
  Dir.mkdir(quant_dir) unless Dir.exist?(quant_dir)

  mask.writeToFile(File.expand_path(base + mask_ext, mask_dir))

  File.open(File.expand_path(base + quant_ext, quant_dir), 'w') do |f|

    f.puts(quant_str)

  end

end