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
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
Applies the STAPLE algorithm to the dataset to determine the true hidden
 segmentation as well as scoring the various segmentations.