C++ image::Header类代码示例

 static Image::Header param (const Image::Header& header, size_t nvols) {
   Image::Header ret (header);
   ret.datatype() = DataType::Float32;
   if (nvols) ret.dim(3) = nvols;
   else ret.set_ndim (3);
   return ret;
 }

      Dynamic::~Dynamic()
      {

        INFO ("Dynamic seeeding required " + str (total_samples) + " samples to draw " + str (total_seeds) + " seeds");

#ifdef DYNAMIC_SEED_DEBUGGING
        const double final_mu = mu();

        // Output seeding probabilites at end of execution
        Image::Header H;
        H.info() = info();
        H.datatype() = DataType::Float32;
        Image::Buffer<float> prob_mean_data ("seed_prob_mean.mif", H), prob_sum_data ("seed_prob_sum.mif", H);
        Image::Buffer<float>::voxel_type prob_mean (prob_mean_data), prob_sum (prob_sum_data);
        VoxelAccessor v (accessor);
        Image::Loop loop;
        for (loop.start (v, prob_mean, prob_sum); loop.ok(); loop.next (v, prob_mean, prob_sum)) {
          if (v.value()) {

            float sum = 0.0;
            size_t count = 0;
            for (Fixel_map<Fixel_TD_seed>::ConstIterator i = begin (v); i; ++i) {
              sum += i().get_seed_prob (final_mu);
              ++count;
            }
            prob_mean.value() = sum / float(count);
            prob_sum .value() = sum;

          }
        }
#endif

      }

void run ()
{
  Image::Buffer<float> dir_buf (argument[0]);
  Image::Buffer<float>::voxel_type dir_vox (dir_buf);

  Image::Header header (argument[0]);
  header.dim(3) = header.dim(3)/3;

  Image::Buffer<float> amp_buf (argument[1], header);
  Image::Buffer<float>::voxel_type amp_vox (amp_buf);

  Image::LoopInOrder loop (dir_vox, "converting directions to amplitudes...", 0, 3);

  for (loop.start (dir_vox, amp_vox); loop.ok(); loop.next (dir_vox, amp_vox)) {
    Math::Vector<float> dir (3);
    dir_vox[3] = 0;
    amp_vox[3] = 0;
    while (dir_vox[3] < dir_vox.dim(3)) {
      dir[0] = dir_vox.value(); ++dir_vox[3];
      dir[1] = dir_vox.value(); ++dir_vox[3];
      dir[2] = dir_vox.value(); ++dir_vox[3];

      float amplitude = 0.0;
      if (std::isfinite (dir[0]) && std::isfinite (dir[1]) && std::isfinite (dir[2]))
        amplitude = Math::norm (dir);

      amp_vox.value() = amplitude;
      ++amp_vox[3];
    }
  }
}

void print_comments (const Image::Header& header)
{
  std::string buffer;
  for (std::vector<std::string>::const_iterator i = header.comments().begin(); i != header.comments().end(); ++i)
    buffer += *i + "\n";
  std::cout << buffer;
}

void run () {

  Image::BufferPreload<bool> input_data (argument[0]);
  Image::BufferPreload<bool>::voxel_type input_voxel (input_data);

  const size_t filter_index = argument[1];

  Image::Filter::Base* filter = NULL;
  switch (filter_index) {
    case 0: filter = create_dilate_filter (input_voxel); break;
    case 1: filter = create_erode_filter  (input_voxel); break;
    case 2: filter = create_lcc_filter    (input_voxel); break;
    case 3: filter = create_median_filter (input_voxel); break;
    default: assert (0);
  }

  Image::Header header;
  header.info() = filter->info();
  Image::Stride::set_from_command_line (header);

  Image::Buffer<bool> output_data (argument[2], header);
  Image::Buffer<bool>::voxel_type output_voxel (output_data);

  filter->set_message (std::string("applying ") + std::string(argument[1]) + " filter to image " + std::string(argument[0]) + "... ");

  switch (filter_index) {
    case 0: (*dynamic_cast<Image::Filter::Dilate*>                    (filter)) (input_voxel, output_voxel); break;
    case 1: (*dynamic_cast<Image::Filter::Erode*>                     (filter)) (input_voxel, output_voxel); break;
    case 2: (*dynamic_cast<Image::Filter::LargestConnectedComponent*> (filter)) (input_voxel, output_voxel); break;
    case 3: (*dynamic_cast<Image::Filter::Median*>                    (filter)) (input_voxel, output_voxel); break;
  }

  delete filter;

}

void print_properties (const Image::Header& header)
{
  std::string buffer;
  for (std::map<std::string, std::string>::const_iterator i = header.begin(); i != header.end(); ++i)
    buffer += i->first + ": " + i->second + "\n";
  std::cout << buffer;
}

void print_vox (const Image::Header& header)
{
  std::string buffer;
  for (size_t i = 0; i < header.ndim(); ++i) {
    if (i) buffer += " ";
    buffer += str (header.vox (i));
  }
  std::cout << buffer << "\n";
}

void print_strides (const Image::Header& header)
{
  std::string buffer;
  std::vector<ssize_t> strides (Image::Stride::get (header));
  Image::Stride::symbolise (strides);
  for (size_t i = 0; i < header.ndim(); ++i) {
    if (i) buffer += " ";
    buffer += header.stride (i) ? str (strides[i]) : "?";
  }
  std::cout << buffer << "\n";
}

inline Point<> get_bounds (const Image::Header& header, Point<int> corner, int axis, const Point<>& ref, const Point<>& d1, const Point<>& d2) {
  const Math::Matrix<float>& T (header.transform());
  Point<> pos (T(0,3), T(1,3), T(2,3));
  for (size_t n = 0; n < 3; ++n) 
    if (corner[n]) 
      pos += header.dim(n) * header.vox(n) * Point<> (T(0,n), T(1,n), T(2,n));
  Point<> dir (T(0,axis), T(1,axis), T(2,axis));
  Point<> ret = get_bounds (pos, dir, ref, d1, d2);
  ret[2] /= -header.vox (axis) * header.dim (axis);
  return ret;
}

void run ()
{
  conv_t conversion = NONE;
  Options opt = get_options ("convert");
  if (opt.size()) {
    switch (int(opt[0][0])) {
      case 0: conversion = OLD; break;
      case 1: conversion = NEW; break;
      case 2:
#ifndef USE_NON_ORTHONORMAL_SH_BASIS
        conversion = NEW;
#else
        conversion = OLD;
#endif
        break;
      case 3: conversion = FORCE_OLDTONEW; break;
      case 4: conversion = FORCE_NEWTOOLD; break;
      default: assert (0); break;
    }
  }

  for (std::vector<ParsedArgument>::const_iterator i = argument.begin(); i != argument.end(); ++i) {

    const std::string path = *i;
    Image::Header H (path);
    if (H.ndim() != 4) {
      WARN ("Image \"" + H.name() + "\" is not 4D and therefore cannot be an SH image");
      continue;
    }
    const size_t lmax = Math::SH::LforN (H.dim(3));
    if (!lmax) {
      WARN ("Image \"" + H.name() + "\" does not contain enough volumes to be an SH image");
      continue;
    }
    if (Math::SH::NforL (lmax) != size_t(H.dim(3))) {
      WARN ("Image \"" + H.name() + "\" does not contain a number of volumes appropriate for an SH image");
      continue;
    }
    if (!H.datatype().is_floating_point()) {
      WARN ("Image \"" + H.name() + "\" does not use a floating-point data type and therefore cannot be an SH image");
      continue;
    }

    if (H.datatype().bytes() == 4)
      check_and_update<float>  (H, conversion);
    else
      check_and_update<double> (H, conversion);

  }

};

void run ()
{
  Image::Header header (argument[0]);

  header.datatype() = DataType::Float32;
  header.set_ndim (4);
  header.dim(3) = 3;
  Image::Stride::set (header, Image::Stride::contiguous_along_axis (3));

  Image::Buffer<float> warp_buffer (argument[1], header);

  Image::ThreadedLoop ("initialising warp image...", warp_buffer, 0, 3)
    .run (write_coordinates (warp_buffer), warp_buffer.voxel());
}

        void oversample_header (Image::Header& header, const std::vector<float>& voxel_size)
        {
          INFO ("oversampling header...");

          Math::Matrix<float> transform (header.transform());
          for (size_t j = 0; j != 3; ++j) {
            for (size_t i = 0; i < 3; ++i)
              header.transform()(i,3) += 0.5 * (voxel_size[j] - header.vox(j)) * transform(i,j);
            header.dim(j) = std::ceil(header.dim(j) * header.vox(j) / voxel_size[j]);
            header.vox(j) = voxel_size[j];
          }
        }

void run ()
{
  Image::Header input_header (argument[0]);
  Image::Buffer<bool> input_data (input_header);
  Image::Buffer<bool>::voxel_type input_voxel (input_data);

  Image::Filter::ConnectedComponents connected_filter(input_voxel);
  Image::Header header (input_data);
  header.info() = connected_filter.info();
  Image::Buffer<int> output_data (argument[1], header);
  Image::Buffer<int>::voxel_type output_vox (output_data);

  Options opt = get_options ("axes");
  std::vector<int> axes;
  if (opt.size()) {
    axes = opt[0][0];
    for (size_t d = 0; d < input_data.ndim(); d++)
      connected_filter.set_ignore_dim (d, true);
    for (size_t i = 0; i < axes.size(); i++) {
      if (axes[i] >= static_cast<int> (input_header.ndim()) || axes[i] < 0)
        throw Exception ("axis supplied to option -ignore is out of bounds");
      connected_filter.set_ignore_dim (axes[i], false);
    }
  }

  opt = get_options ("largest");
  if (opt.size())
    connected_filter.set_largest_only (true);

  opt = get_options ("connectivity");
  if (opt.size())
    connected_filter.set_26_connectivity(true);

  connected_filter.set_message ("computing connected components...");
  connected_filter (input_voxel, output_vox);
}

 Segmented_FOD_receiver (const Image::Header& header, const DWI::Directions::Set& directions) :
   H (header),
   H_fixel (header),
   dirs (directions),
   lmax (0)
 {
   // aPSF does not have data for lmax > 10
   lmax = std::min (size_t(10), Math::SH::LforN (header.dim(3)));
   H.set_ndim (3);
   H.DW_scheme().clear();
   H_fixel.set_ndim (3);
   H_fixel.DW_scheme().clear();
   H_fixel.datatype() = DataType::UInt64;
   H_fixel.datatype().set_byte_order_native();
   H_fixel[Image::Sparse::name_key] = str(typeid(FixelMetric).name());
   H_fixel[Image::Sparse::size_key] = str(sizeof(FixelMetric));
 }

        void generate_header (Image::Header& header, const std::string& tck_file_path, const std::vector<float>& voxel_size)
        {

          Tractography::Properties properties;
          Tractography::Reader<float> file (tck_file_path, properties);

          Streamline<float> tck;
          size_t track_counter = 0;

          Point<float> min_values ( INFINITY,  INFINITY,  INFINITY);
          Point<float> max_values (-INFINITY, -INFINITY, -INFINITY);

          {
            ProgressBar progress ("creating new template image...", 0);
            while (file (tck) && track_counter++ < MAX_TRACKS_READ_FOR_HEADER) {
              for (std::vector<Point<float> >::const_iterator i = tck.begin(); i != tck.end(); ++i) {
                min_values[0] = std::min (min_values[0], (*i)[0]);
                max_values[0] = std::max (max_values[0], (*i)[0]);
                min_values[1] = std::min (min_values[1], (*i)[1]);
                max_values[1] = std::max (max_values[1], (*i)[1]);
                min_values[2] = std::min (min_values[2], (*i)[2]);
                max_values[2] = std::max (max_values[2], (*i)[2]);
              }
              ++progress;
            }
          }

          min_values -= Point<float> (3.0*voxel_size[0], 3.0*voxel_size[1], 3.0*voxel_size[2]);
          max_values += Point<float> (3.0*voxel_size[0], 3.0*voxel_size[1], 3.0*voxel_size[2]);

          header.name() = "tckmap image header";
          header.set_ndim (3);

          for (size_t i = 0; i != 3; ++i) {
            header.dim(i) = std::ceil((max_values[i] - min_values[i]) / voxel_size[i]);
            header.vox(i) = voxel_size[i];
            header.stride(i) = i+1;
            //header.set_units (i, Image::Axis::millimeters);
          }

          //header.set_description (0, Image::Axis::left_to_right);
          //header.set_description (1, Image::Axis::posterior_to_anterior);
          //header.set_description (2, Image::Axis::inferior_to_superior);

          Math::Matrix<float>& M (header.transform());
          M.allocate (4,4);
          M.identity();
          M(0,3) = min_values[0];
          M(1,3) = min_values[1];
          M(2,3) = min_values[2];
          file.close();
        }