<commit_before>/*
 * Copyright (c) 2008-2016 the MRtrix3 contributors
 * 
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/
 * 
 * MRtrix is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * 
 * For more details, see www.mrtrix.org
 * 
 */

#include "header.h"
#include "image_io/default.h"
#include "image_io/mosaic.h"
#include "file/dicom/mapper.h"
#include "file/dicom/image.h"
#include "file/dicom/series.h"
#include "file/dicom/study.h"
#include "file/dicom/patient.h"
#include "file/dicom/tree.h"

namespace MR {
  namespace File {
    namespace Dicom {


      std::unique_ptr<MR::ImageIO::Base> dicom_to_mapper (MR::Header& H, std::vector<std::shared_ptr<Series>>& series)
      {
        assert (series.size() > 0);
        std::unique_ptr<MR::ImageIO::Base> io_handler;

        Patient* patient (series[0]->study->patient);
        std::string sbuf = ( patient->name.size() ? patient->name : "unnamed" );
        sbuf += " " + format_ID (patient->ID);
        if (series[0]->modality.size())
          sbuf += std::string (" [") + series[0]->modality + "]";
        if (series[0]->name.size()) 
          sbuf += std::string (" ") + series[0]->name;
        add_line (H.keyval()["comments"], sbuf);
        H.name() = sbuf;

        // build up sorted list of frames:
        std::vector<Frame*> frames;

        // loop over series list:
        for (const auto series_it : series) {
          try {
            series_it->read();
          }
          catch (Exception& E) { 
            E.display();
            throw Exception ("error reading series " + str (series_it->number) + " of DICOM image \"" + H.name() + "\""); 
          }

          std::sort (series_it->begin(), series_it->end(), compare_ptr_contents());

          // loop over images in each series:
          for (auto image_it : *series_it) {
            // if multi-frame, loop over frames in image:
            if (image_it->frames.size()) {
              std::sort (image_it->frames.begin(), image_it->frames.end(), compare_ptr_contents());
              for (auto frame_it : image_it->frames) 
                frames.push_back (frame_it.get());
            }
            // otherwise add image frame:
            else 
              frames.push_back (image_it.get());
          }
        }

        auto dim = Frame::count (frames);

        if (dim[0]*dim[1]*dim[2] < frames.size())
          throw Exception ("dimensions mismatch in DICOM series");

        if (dim[0]*dim[1]*dim[2] > frames.size())
          throw Exception ("missing image frames for DICOM image \"" + H.name() + "\"");

        if (dim[0] > 1) { // switch axes so slice dim is inner-most:
          std::vector<Frame*> list (frames);
          std::vector<Frame*>::iterator it = frames.begin();
          for (size_t k = 0; k < dim[2]; ++k) 
            for (size_t i = 0; i < dim[0]; ++i) 
              for (size_t j = 0; j < dim[1]; ++j) 
                *(it++) = list[i+dim[0]*(j+dim[1]*k)];
        }

        default_type slice_separation = Frame::get_slice_separation (frames, dim[1]);

        if (series[0]->study->name.size()) 
          add_line (H.keyval()["comments"], std::string ("study: " + series[0]->study->name));

        if (patient->DOB.size()) 
          add_line (H.keyval()["comments"], std::string ("DOB: " + format_date (patient->DOB)));

        if (series[0]->date.size()) {
          sbuf = "DOS: " + format_date (series[0]->date);
          if (series[0]->time.size()) 
            sbuf += " " + format_time (series[0]->time);
          add_line (H.keyval()["comments"], sbuf);
        }

        const Image& image (*(*series[0])[0]);

        if (image.pe_axis != 3 && image.pe_sign) {
          std::string pe_symbol;
          switch (image.pe_axis) {
            case 0: pe_symbol = "i"; break;
            case 1: pe_symbol = "j"; break;
            case 2: pe_symbol = "k"; break;
            default: assert (0);
          }
          if (image.pe_sign < 0)
            pe_symbol += '-';
          add_line (H.keyval()["comments"], "PhaseEncodingDirection: " + pe_symbol);
        }

        if (std::isfinite (image.bandwidth_per_pixel_phase_encode) && image.pe_axis != 3) {
          const default_type effective_echo_spacing = 1.0 / (image.bandwidth_per_pixel_phase_encode * image.dim[image.pe_axis]);
          add_line (H.keyval()["comments"], "EffectiveEchoSpacing: " + str (effective_echo_spacing, 3));
        }

        if (std::isfinite (image.echo_time))
          add_line (H.keyval()["comments"], "EchoTime: " + str (0.001 * image.echo_time, 6));

        size_t nchannels = image.frames.size() ? 1 : image.data_size / (image.dim[0] * image.dim[1] * (image.bits_alloc/8));
        if (nchannels > 1) 
          INFO ("data segment is larger than expected from image dimensions - interpreting as multi-channel data");

        H.ndim() = 3 + (dim[0]*dim[2]>1) + (nchannels>1);

        size_t current_axis = 0;

        if (nchannels > 1) {
          H.stride(3) = 1;
          H.size(3) = nchannels;
          ++current_axis;
        }

        H.stride(0) = ++current_axis;
        H.size(0) = image.dim[0];
        H.spacing(0) = image.pixel_size[0];

        H.stride(1) = ++current_axis;
        H.size(1) = image.dim[1];
        H.spacing(1) = image.pixel_size[1];

        H.stride(2) = ++current_axis;
        H.size(2) = dim[1];
        H.spacing(2) = slice_separation;

        if (dim[0]*dim[2] > 1) {
          H.stride(current_axis) = current_axis+1;
          H.size(current_axis) = dim[0]*dim[2];
          ++current_axis;
        }


        if (image.bits_alloc == 8) 
          H.datatype() = DataType::UInt8;
        else if (image.bits_alloc == 16) {
          H.datatype() = DataType::UInt16;
          if (image.is_BE) 
            H.datatype() = DataType::UInt16 | DataType::BigEndian;
          else 
            H.datatype() = DataType::UInt16 | DataType::LittleEndian;
        }
        else throw Exception ("unexpected number of allocated bits per pixel (" + str (image.bits_alloc) 
            + ") in file \"" + H.name() + "\"");

        H.set_intensity_scaling (image.scale_slope, image.scale_intercept);


        // If multi-frame, take the transform information from the sorted frames; the first entry in the
        // vector should be the first slice of the first volume
        {
          transform_type M;
          const Frame* frame (image.frames.size() ? image.frames[0].get() : &static_cast<const Frame&> (image));

          M(0,0) = -frame->orientation_x[0];
          M(1,0) = -frame->orientation_x[1];
          M(2,0) = +frame->orientation_x[2];

          M(0,1) = -frame->orientation_y[0];
          M(1,1) = -frame->orientation_y[1];
          M(2,1) = +frame->orientation_y[2];

          M(0,2) = -frame->orientation_z[0];
          M(1,2) = -frame->orientation_z[1];
          M(2,2) = +frame->orientation_z[2];

          M(0,3) = -frame->position_vector[0];
          M(1,3) = -frame->position_vector[1];
          M(2,3) = +frame->position_vector[2];

          H.transform() = M;
          std::string dw_scheme = Frame::get_DW_scheme (frames, dim[1], M);
          if (dw_scheme.size())
            H.keyval()["dw_scheme"] = dw_scheme;
        }




        for (size_t n = 1; n < frames.size(); ++n) // check consistency of data scaling:
          if (frames[n]->scale_intercept != frames[n-1]->scale_intercept ||
              frames[n]->scale_slope != frames[n-1]->scale_slope)
            throw Exception ("unable to load series due to inconsistent data scaling between DICOM images");


        if (image.images_in_mosaic) {
          INFO ("DICOM image \"" + H.name() + "\" is in mosaic format");
          if (H.size (2) != 1)
            throw Exception ("DICOM mosaic contains multiple slices in image \"" + H.name() + "\"");

          H.size(0) = image.acq_dim[0];
          H.size(1) = image.acq_dim[1];
          H.size(2) = image.images_in_mosaic;

          if (image.dim[0] % image.acq_dim[0] || image.dim[1] % image.acq_dim[1]) {
            WARN ("acquisition matrix [ " + str (image.acq_dim[0]) + " " + str (image.acq_dim[1]) 
                + " ] does not fit into DICOM mosaic [ " + str (image.dim[0]) + " " + str (image.dim[1]) 
                + " ] -  adjusting matrix size to suit");
            H.size(0) = image.dim[0] / size_t (float(image.dim[0]) / float(image.acq_dim[0]));
            H.size(1) = image.dim[1] / size_t (float(image.dim[1]) / float(image.acq_dim[1]));
          }

          float xinc = H.spacing(0) * (image.dim[0] - H.size(0)) / 2.0;
          float yinc = H.spacing(1) * (image.dim[1] - H.size(1)) / 2.0;
          for (size_t i = 0; i < 3; i++) 
            H.transform()(i,3) += xinc * H.transform()(i,0) + yinc * H.transform()(i,1);

          io_handler.reset (new MR::ImageIO::Mosaic (H, image.dim[0], image.dim[1], H.size (0), H.size (1), H.size (2)));
        }
        else 
          io_handler.reset (new MR::ImageIO::Default (H));

        for (size_t n = 0; n < frames.size(); ++n) 
          io_handler->files.push_back (File::Entry (frames[n]->filename, frames[n]->data));

        return io_handler;
      }



    }
  }
}

<commit_msg>DICOM phase encode info changes<commit_after>/*
 * Copyright (c) 2008-2016 the MRtrix3 contributors
 * 
 * This Source Code Form is subject to the terms of the Mozilla Public
 * License, v. 2.0. If a copy of the MPL was not distributed with this
 * file, You can obtain one at http://mozilla.org/MPL/2.0/
 * 
 * MRtrix is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * 
 * For more details, see www.mrtrix.org
 * 
 */

#include "header.h"
#include "image_io/default.h"
#include "image_io/mosaic.h"
#include "file/dicom/mapper.h"
#include "file/dicom/image.h"
#include "file/dicom/series.h"
#include "file/dicom/study.h"
#include "file/dicom/patient.h"
#include "file/dicom/tree.h"

namespace MR {
  namespace File {
    namespace Dicom {


      std::unique_ptr<MR::ImageIO::Base> dicom_to_mapper (MR::Header& H, std::vector<std::shared_ptr<Series>>& series)
      {
        assert (series.size() > 0);
        std::unique_ptr<MR::ImageIO::Base> io_handler;

        Patient* patient (series[0]->study->patient);
        std::string sbuf = ( patient->name.size() ? patient->name : "unnamed" );
        sbuf += " " + format_ID (patient->ID);
        if (series[0]->modality.size())
          sbuf += std::string (" [") + series[0]->modality + "]";
        if (series[0]->name.size()) 
          sbuf += std::string (" ") + series[0]->name;
        add_line (H.keyval()["comments"], sbuf);
        H.name() = sbuf;

        // build up sorted list of frames:
        std::vector<Frame*> frames;

        // loop over series list:
        for (const auto series_it : series) {
          try {
            series_it->read();
          }
          catch (Exception& E) { 
            E.display();
            throw Exception ("error reading series " + str (series_it->number) + " of DICOM image \"" + H.name() + "\""); 
          }

          std::sort (series_it->begin(), series_it->end(), compare_ptr_contents());

          // loop over images in each series:
          for (auto image_it : *series_it) {
            // if multi-frame, loop over frames in image:
            if (image_it->frames.size()) {
              std::sort (image_it->frames.begin(), image_it->frames.end(), compare_ptr_contents());
              for (auto frame_it : image_it->frames) 
                frames.push_back (frame_it.get());
            }
            // otherwise add image frame:
            else 
              frames.push_back (image_it.get());
          }
        }

        auto dim = Frame::count (frames);

        if (dim[0]*dim[1]*dim[2] < frames.size())
          throw Exception ("dimensions mismatch in DICOM series");

        if (dim[0]*dim[1]*dim[2] > frames.size())
          throw Exception ("missing image frames for DICOM image \"" + H.name() + "\"");

        if (dim[0] > 1) { // switch axes so slice dim is inner-most:
          std::vector<Frame*> list (frames);
          std::vector<Frame*>::iterator it = frames.begin();
          for (size_t k = 0; k < dim[2]; ++k) 
            for (size_t i = 0; i < dim[0]; ++i) 
              for (size_t j = 0; j < dim[1]; ++j) 
                *(it++) = list[i+dim[0]*(j+dim[1]*k)];
        }

        default_type slice_separation = Frame::get_slice_separation (frames, dim[1]);

        if (series[0]->study->name.size()) 
          add_line (H.keyval()["comments"], std::string ("study: " + series[0]->study->name));

        if (patient->DOB.size()) 
          add_line (H.keyval()["comments"], std::string ("DOB: " + format_date (patient->DOB)));

        if (series[0]->date.size()) {
          sbuf = "DOS: " + format_date (series[0]->date);
          if (series[0]->time.size()) 
            sbuf += " " + format_time (series[0]->time);
          add_line (H.keyval()["comments"], sbuf);
        }

        const Image& image (*(*series[0])[0]);

        if (image.pe_axis != 3 && image.pe_sign) {
          std::string pe_symbol;
          switch (image.pe_axis) {
            case 0: pe_symbol = "i"; break;
            case 1: pe_symbol = "j"; break;
            case 2: pe_symbol = "k"; break;
            default: assert (0);
          }
          if (image.pe_sign < 0)
            pe_symbol += '-';
          H.keyval()["PhaseEncodingDirection"] = pe_symbol;
        }

        if (std::isfinite (image.bandwidth_per_pixel_phase_encode) && image.pe_axis != 3) {
          const default_type effective_echo_spacing = 1.0 / (image.bandwidth_per_pixel_phase_encode * image.dim[image.pe_axis]);
          H.keyval()["EffectiveEchoSpacing"] = str (effective_echo_spacing, 3);
          H.keyval()["TotalReadoutTime"] = str (effective_echo_spacing * (image.dim[image.pe_axis] - 1), 3);
        }

        if (std::isfinite (image.echo_time))
          H.keyval()["EchoTime"] = str (0.001 * image.echo_time, 6);

        size_t nchannels = image.frames.size() ? 1 : image.data_size / (image.dim[0] * image.dim[1] * (image.bits_alloc/8));
        if (nchannels > 1) 
          INFO ("data segment is larger than expected from image dimensions - interpreting as multi-channel data");

        H.ndim() = 3 + (dim[0]*dim[2]>1) + (nchannels>1);

        size_t current_axis = 0;

        if (nchannels > 1) {
          H.stride(3) = 1;
          H.size(3) = nchannels;
          ++current_axis;
        }

        H.stride(0) = ++current_axis;
        H.size(0) = image.dim[0];
        H.spacing(0) = image.pixel_size[0];

        H.stride(1) = ++current_axis;
        H.size(1) = image.dim[1];
        H.spacing(1) = image.pixel_size[1];

        H.stride(2) = ++current_axis;
        H.size(2) = dim[1];
        H.spacing(2) = slice_separation;

        if (dim[0]*dim[2] > 1) {
          H.stride(current_axis) = current_axis+1;
          H.size(current_axis) = dim[0]*dim[2];
          ++current_axis;
        }


        if (image.bits_alloc == 8) 
          H.datatype() = DataType::UInt8;
        else if (image.bits_alloc == 16) {
          H.datatype() = DataType::UInt16;
          if (image.is_BE) 
            H.datatype() = DataType::UInt16 | DataType::BigEndian;
          else 
            H.datatype() = DataType::UInt16 | DataType::LittleEndian;
        }
        else throw Exception ("unexpected number of allocated bits per pixel (" + str (image.bits_alloc) 
            + ") in file \"" + H.name() + "\"");

        H.set_intensity_scaling (image.scale_slope, image.scale_intercept);


        // If multi-frame, take the transform information from the sorted frames; the first entry in the
        // vector should be the first slice of the first volume
        {
          transform_type M;
          const Frame* frame (image.frames.size() ? image.frames[0].get() : &static_cast<const Frame&> (image));

          M(0,0) = -frame->orientation_x[0];
          M(1,0) = -frame->orientation_x[1];
          M(2,0) = +frame->orientation_x[2];

          M(0,1) = -frame->orientation_y[0];
          M(1,1) = -frame->orientation_y[1];
          M(2,1) = +frame->orientation_y[2];

          M(0,2) = -frame->orientation_z[0];
          M(1,2) = -frame->orientation_z[1];
          M(2,2) = +frame->orientation_z[2];

          M(0,3) = -frame->position_vector[0];
          M(1,3) = -frame->position_vector[1];
          M(2,3) = +frame->position_vector[2];

          H.transform() = M;
          std::string dw_scheme = Frame::get_DW_scheme (frames, dim[1], M);
          if (dw_scheme.size())
            H.keyval()["dw_scheme"] = dw_scheme;
        }




        for (size_t n = 1; n < frames.size(); ++n) // check consistency of data scaling:
          if (frames[n]->scale_intercept != frames[n-1]->scale_intercept ||
              frames[n]->scale_slope != frames[n-1]->scale_slope)
            throw Exception ("unable to load series due to inconsistent data scaling between DICOM images");


        if (image.images_in_mosaic) {
          INFO ("DICOM image \"" + H.name() + "\" is in mosaic format");
          if (H.size (2) != 1)
            throw Exception ("DICOM mosaic contains multiple slices in image \"" + H.name() + "\"");

          H.size(0) = image.acq_dim[0];
          H.size(1) = image.acq_dim[1];
          H.size(2) = image.images_in_mosaic;

          if (image.dim[0] % image.acq_dim[0] || image.dim[1] % image.acq_dim[1]) {
            WARN ("acquisition matrix [ " + str (image.acq_dim[0]) + " " + str (image.acq_dim[1]) 
                + " ] does not fit into DICOM mosaic [ " + str (image.dim[0]) + " " + str (image.dim[1]) 
                + " ] -  adjusting matrix size to suit");
            H.size(0) = image.dim[0] / size_t (float(image.dim[0]) / float(image.acq_dim[0]));
            H.size(1) = image.dim[1] / size_t (float(image.dim[1]) / float(image.acq_dim[1]));
          }

          float xinc = H.spacing(0) * (image.dim[0] - H.size(0)) / 2.0;
          float yinc = H.spacing(1) * (image.dim[1] - H.size(1)) / 2.0;
          for (size_t i = 0; i < 3; i++) 
            H.transform()(i,3) += xinc * H.transform()(i,0) + yinc * H.transform()(i,1);

          io_handler.reset (new MR::ImageIO::Mosaic (H, image.dim[0], image.dim[1], H.size (0), H.size (1), H.size (2)));
        }
        else 
          io_handler.reset (new MR::ImageIO::Default (H));

        for (size_t n = 0; n < frames.size(); ++n) 
          io_handler->files.push_back (File::Entry (frames[n]->filename, frames[n]->data));

        return io_handler;
      }



    }
  }
}


{'fraction_non_alphanumeric': 0.16708887481528392, 'fraction_numerical': 0.020054887059320244, 'mean_word_length': 2.1565883724804267, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 127, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7287821', 'n_tokens_mistral': 6570, 'n_tokens_neox': 6217, 'n_words': 1404}