00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #ifndef MLN_IO_DICOM_LOAD_HH
00028 # define MLN_IO_DICOM_LOAD_HH
00029
00032
00033 # include <mln/core/image/image2d.hh>
00034 # include <mln/core/image/image3d.hh>
00035
00036 # include <mln/algebra/vec.hh>
00037
00038 # include <gdcm-2.0/gdcmReader.h>
00039 # include <gdcm-2.0/gdcmImageReader.h>
00040 # include <gdcm-2.0/gdcmWriter.h>
00041 # include <gdcm-2.0/gdcmDataSet.h>
00042 # include <gdcm-2.0/gdcmAttribute.h>
00043
00044
00045 namespace mln
00046 {
00047
00048 namespace io
00049 {
00050
00051 namespace dicom
00052 {
00053
00066 template <typename I>
00067 void load(Image<I>& ima,
00068 const std::string& filename);
00069
00070
00071 # ifndef MLN_INCLUDE_ONLY
00072
00073 template <typename V>
00074 inline
00075 image2d<V> load(const std::string& filename)
00076 {
00077 trace::entering("mln::io::gdcm::load");
00078 image2d<V> ima;
00079 trace::exiting("mln::io::gdcm::load");
00080 return ima;
00081 }
00082
00083 template <typename V>
00084 inline
00085 image3d<V> load(const std::string& filename)
00086 {
00087 trace::entering("mln::io::gdcm::load");
00088 image2d<V> ima;
00089 trace::exiting("mln::io::gdcm::load");
00090 return ima;
00091 }
00092
00093
00094 template <typename I>
00095 inline
00096 void load(Image<I>& ima_,
00097 const std::string& filename)
00098 {
00099 trace::entering("mln::io::dicom::load");
00100
00101 I& ima = exact(ima_);
00102
00103 gdcm::ImageReader r;
00104 r.SetFileName(filename.c_str());
00105 if (!r.Read())
00106 {
00107 std::cerr << "error: cannot open file '" << filename << "'!";
00108 abort();
00109 }
00110
00111
00112
00113
00114 gdcm::Image& image = r.GetImage();
00115
00116 char* dataBuffer = new char[image.GetBufferLength()];
00117 image.GetBuffer(dataBuffer);
00118
00119 int ndims = image.GetNumberOfDimensions();
00120 const unsigned int* dims = image.GetDimensions();
00121
00122 unsigned short bits_allocated = image.GetPixelFormat().GetBitsAllocated();
00123 unsigned short bytes_allocated = bits_allocated / 8;
00124 unsigned short bits_stored = image.GetPixelFormat().GetBitsStored();
00125 unsigned short samples_per_pixel = image.GetPixelFormat().GetSamplesPerPixel();
00126
00127 unsigned int offset = 8 - (bits_allocated - bits_stored);
00128 unsigned int off_pow = 1;
00129 for (unsigned int i = 0; i < offset; ++i)
00130 {
00131 off_pow *= 2;
00132 }
00133
00134 if (mln_site_(I)::dim != ndims)
00135 {
00136 std::cerr << "error: dimension mismatch" << std::endl;
00137 abort();
00138 }
00139
00140 algebra::vec<mln_site_(I)::dim, unsigned int> vmin;
00141 algebra::vec<mln_site_(I)::dim, unsigned int> vmax;
00142 algebra::vec<mln_site_(I)::dim, unsigned int> vdims;
00143 int j = ndims - 1;
00144 for (int i = 0; i < ndims; ++i, --j)
00145 {
00146 vmin[i] = 0;
00147 vmax[i] = dims[i] - 1;
00148 if (i == 0)
00149 vdims[j] = 1;
00150 else
00151 vdims[j] = dims[i - 1] * vdims[j + 1];
00152 }
00153
00154 if (ndims > 1)
00155 {
00156 std::swap(vmin[0], vmin[1]);
00157 std::swap(vmax[0], vmax[1]);
00158 }
00159
00160 mln_site(I) pmin(vmin);
00161 mln_site(I) pmax(vmax);
00162 mln_concrete(I) result(box<mln_site(I)>(pmin, pmax));
00163 initialize(ima, result);
00164 mln_piter(I) p(ima.domain());
00165 unsigned int index = 0;
00166
00167 for_all(p)
00168 {
00169 index = 0;
00170 for (int i = 0; i < ndims; ++i)
00171 {
00172 index += p[i] * vdims[i];
00173 }
00174
00175 mln_value(I) v = (unsigned char) dataBuffer[(index * bytes_allocated) * samples_per_pixel];
00176
00177 for (unsigned int j = 0; j < bytes_allocated; ++j)
00178 {
00179 v += ((unsigned char) dataBuffer[(index * bytes_allocated + j) * samples_per_pixel]) * 256 * j;
00180 }
00181
00182 ima(p) = v;
00183 }
00184
00185 delete(dataBuffer);
00186
00187 trace::exiting("mln::io::dicom::load");
00188 }
00189
00190 # endif // ! MLN_INCLUDE_ONLY
00191
00192 }
00193
00194 }
00195
00196 }
00197
00198
00199 #endif // ! MLN_IO_DICOM_LOAD_HH