Class: RTKIT::Staple

Inherits:
Object
  • Object
show all
Defined in:
lib/rtkit/staple.rb

Overview

The Staple class is used for simultaneously evaluating the performance of multiple volume segmentations (typically derived from a RT Structure Set) as well as establishing the hidden true segmentation based on probabilistic analysis of the supplied rater decisions. The determined true segmentation can easily be exported to a RT Structure Set for external use.

THEORY

Complete data: (D, T) Probability mass function of the complete data: f(D,T | p,q) Task: Which performance level parameters (p,q) will maximize the complete data log likelihood function: (p’,q’) = arg max_pq ln f(D,T | p,q)

Indices: D Voxel nr: i Segmentation nr: j Iteration nr: k

The Expectation-Maximization algorithm approaches the problem of maximizing the incomplete data log likelihood equation by proceeding iteratively with estimation and maximization of the complete data log likelihood function.

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(bin_matcher, options = {}) ⇒ Staple

Creates a Staple instance for the provided segmented volumes.

Parameters

  • bin_matcher – An BinMatcher instance containing at least two volumes.

  • options – A hash of parameters.

Options

  • :max_iterations – Integer. The maximum number of iterations to use in the STAPLE algorithm. Defaults to 25.

Raises:

  • (ArgumentError) 


60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
 
# File 'lib/rtkit/staple.rb', line 60

def initialize(bin_matcher, options={})
  raise ArgumentError, "Invalid argument 'bin_matcher'. Expected BinMatcher, got #{bin_matcher.class}." unless bin_matcher.is_a?(BinMatcher)
  raise ArgumentError, "Invalid argument 'bin_matcher'. Expected BinMatcher with at least 2 volumes, got #{bin_matcher.volumes.length}." unless bin_matcher.volumes.length > 1
  # Verify that the volumes have equal dimensions (columns and rows):
  volumes = bin_matcher.narrays(sort=false)
  raise ArgumentError, "Invalid argument 'bin_matcher'. Expected BinMatcher with volumes having the same number of columns, got #{volumes.collect{|v| v.shape[1]}.uniq}." unless volumes.collect{|v| v.shape[1]}.uniq.length == 1
  raise ArgumentError, "Invalid argument 'bin_matcher'. Expected BinMatcher with volumes having the same number of rows, got #{volumes.collect{|v| v.shape[2]}.uniq}." unless volumes.collect{|v| v.shape[2]}.uniq.length == 1
  # Make sure the volumes of the BinMatcher instance are comparable:
  bin_matcher.fill_blanks
  bin_matcher.sort_volumes
  @volumes = bin_matcher.narrays(sort=false)
  @original_volumes = @volumes.dup
  # Verify that the volumes have the same number of frames:
  raise ArgumentError, "Invalid argument 'bin_matcher'. Expected BinMatcher with volumes having the same number of frames, got #{@volumes.collect{|v| v.shape[0]}.uniq}." unless @volumes.collect{|v| v.shape[0]}.uniq.length == 1
  @bm = bin_matcher
  # Options:
  @max_iterations = options[:max_iterations] || 25
end

Instance Attribute Details

#decisions ⇒ Object (readonly)

An NArray containing all rater decisions (dimensions n*r).



29
30
31
 
# File 'lib/rtkit/staple.rb', line 29

def decisions
  @decisions
end

#max_iterations ⇒ Object

The maximum number of iterations to use in the STAPLE algorithm.



31
32
33
 
# File 'lib/rtkit/staple.rb', line 31

def max_iterations
  @max_iterations
end

#n ⇒ Object (readonly)

Number of voxels in the volume to be evaluated.



33
34
35
 
# File 'lib/rtkit/staple.rb', line 33

def n
  @n
end

#p ⇒ Object (readonly)

Sensitivity float vector (length r). Each index contains a score from 0 (worst) to 1 (all true voxels segmented by the rater).



35
36
37
 
# File 'lib/rtkit/staple.rb', line 35

def p
  @p
end

#phi ⇒ Object (readonly)

An NArray containing the results of the Staple analysis (dimensions 2*r).



37
38
39
 
# File 'lib/rtkit/staple.rb', line 37

def phi
  @phi
end

#q ⇒ Object (readonly)

Specificity float vector (length r). Each index contains a score from 0 (worst) to 1 (none of the true remaining voxels are segmented by the rater).



39
40
41
 
# File 'lib/rtkit/staple.rb', line 39

def q
  @q
end

#r ⇒ Object (readonly)

Number of raters to be evaluated.



41
42
43
 
# File 'lib/rtkit/staple.rb', line 41

def r
  @r
end

#true_segmentation ⇒ Object (readonly)

An NArray containing the determined true segmentation (dimensions equal to that of the input volumes).



43
44
45
 
# File 'lib/rtkit/staple.rb', line 43

def true_segmentation
  @true_segmentation
end

#vectors ⇒ Object (readonly)

The decision vectors used (derived from the supplied volumes).



45
46
47
 
# File 'lib/rtkit/staple.rb', line 45

def vectors
  @vectors
end

#weights ⇒ Object (readonly)

A float vector containing the weights assigned to each voxel (when rounded becomes the true segmentation) (length n).



47
48
49
 
# File 'lib/rtkit/staple.rb', line 47

def weights
  @weights
end

Instance Method Details

#==(other) ⇒ Object Also known as: eql?

Returns true if the argument is an instance with attributes equal to self.



81
82
83
84
85
 
# File 'lib/rtkit/staple.rb', line 81

def ==(other)
  if other.respond_to?(:to_staple)
    other.send(:state) == state
  end
end

#hash ⇒ Object

Generates a Fixnum hash value for this instance.



91
92
93
 
# File 'lib/rtkit/staple.rb', line 91

def hash
  state.hash
end

#remove_empty_indices ⇒ Object

Along each dimension of the input volume, removes any index (slice, column or row) which is empty in all volumes. The result is a reduced volume used for the analysis, yielding scores with better contrast on specificity. This implementation aims to be independent of the number of dimensions in the input segmentation.



99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
 
# File 'lib/rtkit/staple.rb', line 99

def remove_empty_indices
  # It only makes sense to run this volume reduction if the number of dimensions are 2 or more:
  if @original_volumes.first.dim > 1
    # To be able to reconstruct the volume later on, we need to keep track of the original indices
    # of the indices that remain in the new, reduced volume:
    @original_indices = Array.new(@original_volumes.first.dim)
    # For a volume the typical meaning of the dimensions will be: slice, column, row
    @original_volumes.first.shape.each_with_index do |size, dim_index|
      extract = Array.new(@original_volumes.first.dim, true)
      segmented_indices = Array.new
      size.times do |i|
        extract[dim_index] = i
        segmented = false
        @volumes.each do |volume|
          segmented = true if volume[*extract].max > 0
        end
        segmented_indices << i if segmented
      end
      @original_indices[dim_index] = segmented_indices
      extract[dim_index] = segmented_indices
      # Iterate each volume and pull out segmented indices:
      if segmented_indices.length < size
        @volumes.collect!{|volume| volume = volume[*extract]}
      end
    end
  end
end

#solve ⇒ Object

Applies the STAPLE algorithm to the dataset to determine the true hidden segmentation as well as scoring the various segmentations.



130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
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
202
203
204
 
# File 'lib/rtkit/staple.rb', line 130

def solve
  set_parameters
  # Vectors holding the values used for calculating the weights:
  a = NArray.float(@n)
  b = NArray.float(@n)
  # Set an initial estimate for the probabilities of true segmentation:
  @n.times do |i|
    @weights_current[i] = @decisions[i, true].mean
  end
  # Proceed iteratively until we have converged to a local optimum:
  k = 0
  while k < max_iterations do
    # Copy weights:
    @weights_previous = @weights_current.dup
    # E-step: Estimation of the conditional expectation of the complete data log likelihood function.
    # Deriving the estimator for the unobserved true segmentation (T).
    @n.times do |i|
      voxel_decisions = @decisions[i, true]
      # Find the rater-indices for this voxel where the raters' decisions equals 1 and 0:
      positive_indices, negative_indices = (voxel_decisions.eq 1).where2
      # Determine ai:
      # Multiply by corresponding sensitivity (or 1 - sensitivity):
      a_decision1_factor = (positive_indices.length == 0 ? 1 : @p[positive_indices].prod)
      a_decision0_factor = (negative_indices.length == 0 ? 1 : (1 - @p[negative_indices]).prod)
      a[i] = @weights_previous[i] * a_decision1_factor * a_decision0_factor
      # Determine bi:
      # Multiply by corresponding specificity (or 1 - specificity):
      b_decision0_factor = (negative_indices.length == 0 ? 1 : @q[negative_indices].prod)
      b_decision1_factor = (positive_indices.length == 0 ? 1 : (1 - @q[positive_indices]).prod)
      b[i] = @weights_previous[i] * b_decision0_factor * b_decision1_factor
      # Determine Wi: (take care not to divide by zero)
      if a[i] > 0 or b[i] > 0
        @weights_current[i] = a[i] / (a[i] + b[i])
      else
        @weights_current[i] = 0
      end
    end
    # M-step: Estimation of the performance parameters by maximization.
    # Finding the values of the expert performance level parameters that maximize the conditional expectation
    # of the complete data log likelihood function (phi - p,q).
    @r.times do |j|
      voxel_decisions = @decisions[true, j]
      # Find the voxel-indices for this rater where the rater's decisions equals 1 and 0:
      positive_indices, negative_indices = (voxel_decisions.eq 1).where2
      # Determine sensitivity:
      # Sum the weights for the indices where the rater's decision equals 1:
      sum_positive = (positive_indices.length == 0 ? 0 : @weights_current[positive_indices].sum)
      @p[j] = sum_positive / @weights_current.sum
      # Determine specificity:
      # Sum the weights for the indices where the rater's decision equals 0:
      sum_negative = (negative_indices.length == 0 ? 0 : (1 - @weights_current[negative_indices]).sum)
      @q[j] = sum_negative / (1 - @weights_current).sum
    end
    # Bump our iteration index:
    k += 1
    # Abort if we have reached the local optimum: (there is no change in the sum of weights)
    if @weights_current.sum - @weights_previous.sum == 0
#puts "Iteration aborted as optimum solution was found!" if @verbose
#logger.info("Iteration aborted as optimum solution was found!")
      break
    end
  end
  # Set the true segmentation:
  @true_segmentation_vector = @weights_current.round
  # Set the weights attribute:
  @weights = @weights_current
  # As this vector doesn't make much sense to the user, it must be converted to a volume. If volume reduction has
  # previously been performed, this must be taken into account when transforming it to a volume:
  construct_segmentation_volume
  # Construct a BinVolume instance for the true segmentation and add it as a master volume to the BinMatcher instance.
  update_bin_matcher
  # Set the phi variable:
  @phi[0, true] = @p
  @phi[1, true] = @q
end

#to_staple ⇒ Object

Returns self.



208
209
210
 
# File 'lib/rtkit/staple.rb', line 208

def to_staple
  self
end