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 collapse

DEFAULT_SEG_CH = 
2
DEFAULT_SEG_PL = 
8
DEFAULT_BEAD_RADIUS = 
24.0
DEFAULT_THREADS = 
1
VERSION = 
'0.9.6'

Class Method Summary collapse

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.



209
210
211
 
# File 'lib/find_beads.rb', line 209

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



222
223
224
 
# File 'lib/find_beads.rb', line 222

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.



65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
 
# File 'lib/find_beads.rb', line 65

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

.compute_correlation(mask, im, ch0, ch1, cens, opts) ⇒ Hash

Runs a correlation analysis between two channels of the image. This is done by directly computing the correlation for each bead and then computing an effective background by swapping the x- and y- axes of one channel within each bead.

Parameters:

  • mask (Image) —

    the mask of the beads

  • ch0 (Integer) —

    the first channel for the correlation

  • ch1 (Integer) —

    the second channel for the correlation

  • cens (Hash) —

    a hash mapping bead labels in the mask to the centroids of the beads

  • opts (Hash) —

    the options hash

Returns:

  • (Hash) —

    a hash mapping :norm_corr, :corr, and :bg_corr to a hash mapping bead labels to their measurements



427
428
429
430
431
432
433
434
435
436
 
# File 'lib/find_beads.rb', line 427

def self.compute_correlation(mask, im, ch0, ch1, cens, opts)
  result = {norm_corr: {}, corr: {}, bg_corr: {}}
  cens.each do |id, cen|
    corr_for_bead = compute_single_bead_correlation(mask, im, ch0, ch1, cen, id, opts[:beadradius], true)
    result[:norm_corr][id] = corr_for_bead[:norm_corr]
    result[:corr][id] = corr_for_bead[:corr]
    result[:bg_corr][id] = corr_for_bead[:bg_corr]
  end
  result
end

.compute_single_bead_correlation(mask, im, ch0, ch1, cen, id, bead_radius, do_normalization) ⇒ Hash

Computes the two-channel correlation for a single bead.

Parameters:

  • mask (Image) —

    the mask of labeled, segmented beads

  • im (Image) —

    the image of the beads (single z-section, multiple channels)

  • ch0 (Integer) —

    the (0-based) index of the first channel of the correlation

  • ch1 (Integer) —

    the (0-based) index of the second channel of the correlation

  • cen (Array) —

    an array containing the x,y coordinates of the centroid of the current bead

  • id (Integer) —

    the label of the current bead in the mask

  • bead_radius (Numeric) —

    the radius of a bead

  • do_normalization (Boolean) —

    whether to normalize to the geometric mean of the autocorrelations of the two channels

Returns:

  • (Hash) —

    a hash containing components :norm_corr (the background-subtracted, normalized correlation), :corr (the non-normalized, non-subtracted two-channel correlation), and :bg_corr (the non-normalized background two-channel correlation, calculated by rotating one channel of the bead 90 degrees around the bead centroid)



371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
 
# File 'lib/find_beads.rb', line 371

def self.compute_single_bead_correlation(mask, im, ch0, ch1, cen, id, bead_radius, do_normalization)
  normalization = if do_normalization then
                    auto_00 = compute_single_bead_correlation(mask, im, ch0, ch0, cen, id, bead_radius, false)[:norm_corr]
                    auto_11 = compute_single_bead_correlation(mask, im, ch1, ch1, cen, id, bead_radius, false)[:norm_corr]
                    Math.sqrt(auto_00 * auto_11)*(auto_00.abs/auto_00)
                  else
                    1
                  end

  box_lower = ImageCoordinate[cen[0] - bead_radius - 1, cen[1] - bead_radius - 1,0,0,0]
  box_upper = ImageCoordinate[cen[0] + bead_radius + 2, cen[1] + bead_radius + 2,1,1,1]
  mask.setBoxOfInterest(box_lower, box_upper)
  
  ch_coord = ImageCoordinate[0,0,0,0,0]
  corr_sum = 0
  bg_sum = 0
  count = 0
  bg_count = 0

  mask.each do |ic|
    next unless mask[ic] == id
    ch_coord.setCoord(ic)
    ch_coord[:c] = ch0
    val = im[ch_coord]
    ch_coord[:c] = ch1
    corr_sum += val * im[ch_coord]
    count += 1
    mirror_coord!(ch_coord, cen[0], cen[1])
    if im.inBounds(ch_coord) then
      bg_sum += val * im[ch_coord]
      bg_count += 1
    end
  end
    
  mask.clearBoxOfInterest
  box_lower.recycle
  box_upper.recycle
  ch_coord.recycle

  {norm_corr: (corr_sum/count - bg_sum/bg_count)/normalization , corr: corr_sum/count, bg_corr: bg_sum/bg_count}

end

.format_correlation_output(corr_output) ⇒ String

Organizes the correlation data into a csv string where each row is a bead and each column is a measurement. Also creates a header row.

Parameters:

  • corr_output (Hash) —

    a hash formatted like the output of FindBeads.compute_correlation

Returns:

  • (String) —

    a string containing the same data in csv format.



475
476
477
478
479
480
481
482
483
484
485
 
# File 'lib/find_beads.rb', line 475

def self.format_correlation_output(corr_output)
  header_row = ["bead index", "normalized, corrected correlation", "correlation", "background correlation"]
  
  CSV.generate do |csv|
    csv << header_row
    
    corr_output[:corr].each_key do |id|
      csv << [id, corr_output[:norm_corr][id], corr_output[:corr][id], corr_output[:bg_corr][id]]
    end
  end
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.



118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
 
# File 'lib/find_beads.rb', line 118

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)

  (next_dist_proj - closest_dist_proj < cutoff)
end

.mask_from_image(im, opts, centroid_storage = nil) ⇒ Object

Generates a segmented mask of beads from an image.

Parameters:

  • im (Image) —

    the image to segment

  • opts (Hash) —

    a hash of commandline arguments.

  • centroid_storage (OpenStruct, #centroids=) (defaults to: nil) —

    an optional object to which the bead centroids will be fed using #centroids= default=nil, causes them not to be set



235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
 
# File 'lib/find_beads.rb', line 235

def self.mask_from_image(im, opts, centroid_storage = nil)
  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)

  remapped_cens = {}
  cen_coord = ImageCoordinate[0,0,0,0,0]

  cens.each do |k, cen|
    cen_coord[:x] = cen[0]
    cen_coord[:y] = cen[1]
    val = final_mask[cen_coord]
    
    if val > 0 then
      remapped_cens[val] = cen
    end
  end

  cen_coord.recycle

  if centroid_storage then
    centroid_storage.centroids= remapped_cens
  end

  final_mask

end

.mirror_coord!(coord, cen_x, cen_y) ⇒ ImageCoordinate

Swaps the x and y coordinates of an ImageCoordinate in place relative to a supplied origin.

Parameters:

  • coord (ImageCoordinate) —

    the coordinate whose components will be swapped.

  • cen_x (Numeric) —

    the x-coordinate of the origin relative to which the swap will be calculated

  • cen_y (Numeric) —

    the y-coordinate of the origin relative to which the swap will be calculated

Returns:

  • (ImageCoordinate) —

    coord



347
348
349
350
351
352
353
 
# File 'lib/find_beads.rb', line 347

def self.mirror_coord!(coord, cen_x, cen_y)
  rel_x = coord[:x] - cen_x
  rel_y = coord[:y] - cen_y
  coord[:x] = (rel_y + cen_x).round.to_i
  coord[:y] = (rel_x + cen_y).round.to_i
  coord
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.



493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
 
# File 'lib/find_beads.rb', line 493

def self.process_file(fn, opts=nil)
  puts "processing #{fn}"

  cens = if opts and opts[:correlation_channels] then
           OpenStruct.new
         else
           nil
         end

  im = RImageAnalysisTools.get_image(fn)
  mask = mask_from_image(im, opts, cens)
  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

  outdat = if opts and opts[:correlation_channels] then
             q = compute_correlation(mask, proj, opts[:correlation_channels][0], opts[:correlation_channels][1], cens.centroids, opts)
             format_correlation_output(q)
           else
             met = IntensityPerPixelMetric.new
             q = met.quantify(mask, is)
             Java::edu.stanford.cfuller.imageanalysistools.frontend.LocalAnalysis.generateDataOutputString(q, nil)
           end

  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)



100
101
102
103
104
105
 
# File 'lib/find_beads.rb', line 100

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.



160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
 
# File 'lib/find_beads.rb', line 160

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.



527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
 
# File 'lib/find_beads.rb', line 527

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
    opt :correlation_channels, "Runs correlation between the specified two comma separated channels", :type => :string, :default => nil
  end

  if opts[:correlation_channels] then
    opts[:correlation_channels] = opts[:correlation_channels].gsub("\s+", "").split(",").map(&:to_i)
  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



446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
 
# File 'lib/find_beads.rb', line 446

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