15 #include "gdcmImageReader.h"
20 for(
int i=0; i<s.size(); i++){
21 result[i]=tolower(s[i]);
28 inline Grid3<T>::Grid3(
int width,
int height,
int depth,
double x_spacing,
double y_spacing,
double z_spacing){
36 data=std::vector<T>(NumberOfElements());
47 data=std::vector<T>(NumberOfElements());
51 inline void Grid3<T>::Init(
int width,
int height,
int depth,
double x_spacing,
double y_spacing,
double z_spacing) {
58 data=std::vector<T>(NumberOfElements());
79 result.
x = p.
x / GetSpacingX();
80 result.
y = p.
y / GetSpacingY();
81 result.
z = p.
z / GetSpacingZ();
89 result.
x = p.
x * GetSpacingX();
90 result.
y = p.
y * GetSpacingY();
91 result.
z = p.
z * GetSpacingZ();
97 return _width*_height*_depth;
137 return std::string(
typeid(T).name() );
143 return x+ y*_width + z*_width*_height;
149 c.
z = index/(_width*_height);
150 index = index%(_width*_height);
158 return (x>=0 && x<_width && y>=0 && y<_height && z>=0 && z<_depth);
163 return IsInside(p.
x,p.
y,p.
z);
168 return background_value;
187 else if(mode==
CONSTANT && !IsInside(x,y,z)){
188 return background_value;
191 return data[Offset(x,y,z)];
199 p.
x=min(p.
x,_width-1);
201 p.
y=min(p.
y,_height-1);
203 p.
z=min(p.
z,_depth-1);
205 else if(mode==
CONSTANT && !IsInside(p)){
206 return background_value;
209 return data[Offset(p.
x,p.
y,p.
z)];
214 double xt=x-floor(x);
215 double yt=y-floor(y);
216 double zt=z-floor(z);
218 int x1=int(floor(x));
220 int y1=int(floor(y));
222 int z1=int(floor(z));
225 return (1-zt)*((1-yt)*((1-xt)*operator()(x1,y1,z1,
true)+
226 (xt)*
operator()(x2,y1,z1,
true,mode))+
227 (yt)*((1-xt)*operator()(x1,y2,z1,
true,mode)+
228 (xt)*
operator()(x2,y2,z1,
true,mode)))+
229 (zt)*((1-yt)*((1-xt)*operator()(x1,y1,z2,
true,mode)+
230 (xt)*
operator()(x2,y1,z2,
true,mode))+
231 (yt)*((1-xt)*operator()(x1,y2,z2,
true,mode)+
232 (xt)*
operator()(x2,y2,z2,
true,mode)));
237 return LinearAt(p.
x,p.
y,p.
z, mode);
244 Vec3i coord=Coordinate(index);
245 if(coord.
x+1<_width){
246 n.push_back(Offset(coord.
x+1, coord.
y, coord.
z));
250 n.push_back(Offset(coord.
x-1, coord.
y, coord.
z));
253 if(coord.
y+1<_height){
254 n.push_back(Offset(coord.
x, coord.
y+1, coord.
z));
258 n.push_back(Offset(coord.
x, coord.
y-1, coord.
z));
261 if(coord.
z+1<_depth){
262 n.push_back(Offset(coord.
x, coord.
y, coord.
z+1));
266 n.push_back(Offset(coord.
x, coord.
y, coord.
z-1));
276 n.push_back(
Vec3i(p.
x+1, p.
y, p.
z));
280 n.push_back(
Vec3i(p.
x-1, p.
y, p.
z));
284 n.push_back(
Vec3i(p.
x, p.
y+1, p.
z));
288 n.push_back(
Vec3i(p.
x, p.
y-1, p.
z));
292 n.push_back(
Vec3i(p.
x, p.
y, p.
z+1));
296 n.push_back(
Vec3i(p.
x, p.
y, p.
z-1));
302 std::fill(data.begin(),data.end(),value);
307 return *std::max_element(data.begin(), data.end());
312 return *std::min_element(data.begin(), data.end());
361 std::copy(src.
data.begin(), src.
data.end(), data.begin());
368 std::swap_ranges(v.
data.begin(), v.
data.end(), data.begin());
374 if(filename.find(
".vtk")!=std::string::npos){
375 return ReadVTK(filename);
377 else if(filename.find(
".mhd")!=string::npos){
378 return ReadMetaImage(filename);
380 else if(filename.find(
".mha")!=string::npos){
381 return ReadMetaImage(filename);
383 else if(filename.find(
".nii")!=string::npos){
384 return ReadNIFTI(filename);
386 else if(filename.find(
".dcm")!=string::npos){
387 return ReadDICOM(filename);
411 is.open( filename.c_str(), ios::in | ios::binary );
418 std::stringstream ss;
425 std::string checkvtk;
428 if( !(
ToLower(checkvtk) ==
"vtk" )){
return false; }
433 bool read_dimensions=
false;
435 bool read_type=
false;
436 bool read_asciibinary=
false;
437 bool read_origin=
false;
443 double spacing_x=1,spacing_y=1,spacing_z=1;
444 double origin_x=0,origin_y=0,origin_z=0;
447 while(!done && getline(is,line)){
449 std::stringstream ss2(line);
452 if(
ToLower(keyword)==
"dimensions"){
457 read_dimensions=
true;
459 else if(
ToLower(keyword)==
"binary"){
460 read_asciibinary=
true;
462 else if(
ToLower(keyword)==
"ascii"){
464 read_asciibinary=
true;
466 else if(
ToLower(keyword)==
"dataset"){
468 if(!(
ToLower(dataset_type)==
"structured_points")){
472 else if(
ToLower(keyword)==
"spacing" ||
ToLower(keyword)==
"aspect_ratio" ){
477 else if(
ToLower(keyword)==
"origin"){
483 else if(
ToLower(keyword)==
"scalars"){
488 else if(
ToLower(keyword)==
"lookup_table"){
492 done=read_dimensions && read_lut && read_type && read_asciibinary;
509 metadata.SetDouble(
"ORIGIN_X",origin_x);
510 metadata.SetDouble(
"ORIGIN_Y",origin_y);
511 metadata.SetDouble(
"ORIGIN_Z",origin_z);
523 if(typestring==
"unsigned_char"){
524 uint8_t* tmpdata=
new uint8_t[nelem];
525 is.read((
char*)tmpdata, nelem*
sizeof(uint8_t));
526 #pragma omp parallel for
527 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
529 metadata.SetString(
"ORIGINAL_TYPE",
"UINT8");
532 else if(typestring==
"unsigned_short"){
533 uint16_t* tmpdata=
new uint16_t[nelem];
535 #pragma omp parallel for
536 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
538 metadata.SetString(
"ORIGINAL_TYPE",
"UINT16");
541 else if(typestring==
"unsigned_int"){
542 uint32_t* tmpdata=
new uint32_t[nelem];
544 #pragma omp parallel for
545 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
547 metadata.SetString(
"ORIGINAL_TYPE",
"UINT32");
550 else if(typestring==
"char"){
551 int8_t* tmpdata=
new int8_t[nelem];
553 #pragma omp parallel for
554 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
556 metadata.SetString(
"ORIGINAL_TYPE",
"INT8");
559 else if(typestring==
"short"){
560 int16_t* tmpdata=
new int16_t[nelem];
562 #pragma omp parallel for
563 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
565 metadata.SetString(
"ORIGINAL_TYPE",
"INT16");
568 else if(typestring==
"int"){
569 int32_t* tmpdata=
new int32_t[nelem];
571 #pragma omp parallel for
572 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
574 metadata.SetString(
"ORIGINAL_TYPE",
"INT32");
577 else if(typestring==
"float"){
578 float* tmpdata=
new float[nelem];
580 #pragma omp parallel for
581 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
583 metadata.SetString(
"ORIGINAL_TYPE",
"FLOAT");
586 else if(typestring==
"double"){
587 double* tmpdata=
new double[nelem];
589 #pragma omp parallel for
590 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
592 metadata.SetString(
"ORIGINAL_TYPE",
"DOUBLE");
617 outfile.open(filename,ofstream::out|ofstream::binary);
618 if(!outfile) {
return false; }
620 long nelem= _width*_height*_depth;
622 double origin_x=metadata.GetDouble(
"ORIGIN_X",0);
623 double origin_y=metadata.GetDouble(
"ORIGIN_Y",0);
624 double origin_z=metadata.GetDouble(
"ORIGIN_Z",0);
626 outfile<<
"# vtk DataFile Version 3.0"<<endl;
627 outfile<<
"Created using the Grid3 library"<<endl;
629 string Datamode=
"BINARY";
631 outfile<<Datamode<<endl;
633 outfile<<
"DATASET STRUCTURED_POINTS"<<endl;
635 outfile<<
"DIMENSIONS "<<_width<<
" "<<_height<<
" "<<_depth<<endl;
636 outfile<<
"ORIGIN "<<origin_x<<
" "<<origin_y<<
" "<<origin_z<<endl;
637 outfile<<
"SPACING "<<_spacingX<<
" "<<_spacingY<<
" "<<_spacingZ<<endl;
638 outfile<<
"POINT_DATA "<<nelem<<endl;
640 string datatype=GetDataType();
641 string vtk_datatype=
"";
643 if(datatype==
"UINT8"){vtk_datatype=
"unsigned_char";}
644 else if(datatype==
"INT8"){vtk_datatype=
"char";}
645 else if(datatype==
"UINT16"){vtk_datatype=
"unsigned_short";}
646 else if(datatype==
"INT16"){vtk_datatype=
"short";}
647 else if(datatype==
"UINT32"){vtk_datatype=
"unsigned_int";}
648 else if(datatype==
"INT32"){vtk_datatype=
"int";}
649 else if(datatype==
"float"){vtk_datatype=
"float";}
650 else if(datatype==
"double"){vtk_datatype=
"double";}
652 outfile<<
"SCALARS image_data "<<vtk_datatype<<endl;
653 outfile<<
"LOOKUP_TABLE default"<<endl;
681 map<string,string> header_info;
685 if(!f){
return false;}
687 string line,key, value;
692 size_t pos=line.find(
"=");
693 if(pos!=string::npos){
695 key=line.substr(0,pos);
696 value=line.substr(pos);
699 size_t first_non_space=key.find_first_not_of(
" ");
700 size_t last_non_space=key.find_last_not_of(
" ");
701 key=key.substr(first_non_space,last_non_space+1);
704 first_non_space=value.find_first_not_of(
" ");
705 last_non_space=value.find_last_not_of(
" ");
706 value=value.substr(first_non_space,last_non_space+1);
708 header_info[key]=value;
714 map<string,string>::iterator it;
718 it=header_info.find(
"ObjectType");
719 if(it==header_info.end() || it->second!=
"Image"){
return false;}
722 it=header_info.find(
"NDims");
723 if(it==header_info.end()){
return false;}
725 ss=stringstream(it->second);
727 if(!(ndims==2 || ndims==3)){
return false;}
731 it=header_info.find(
"DimSize");
732 if(it==header_info.end()){
return false;}
733 ss=stringstream(it->second);
734 int dimx, dimy, dimz;
739 if(ndims==3){ss>>dimz;}
743 it=header_info.find(
"ElementDataFile");
744 if(it==header_info.end()){
return false;}
745 string datafile=it->second;
748 it=header_info.find(
"ElementType");
749 if(it==header_info.end()){
return false;}
750 string typestring=it->second;
753 bool BinaryData=
true;
754 it=header_info.find(
"BinaryData");
755 if(it!=header_info.end()){
756 BinaryData=(it->second==
"True");
761 it=header_info.find(
"BinaryDataByteOrderMSB");
762 if(it!=header_info.end()){
763 BinaryDataByteOrderMSB=(it->second==
"True");
768 it=header_info.find(
"HeaderSize");
769 if(it!=header_info.end()){
770 ss=stringstream(it->second);
775 it=header_info.find(
"ElementSize");
776 double voxelsize_x=1, voxelsize_y=1, voxelsize_z=1;
777 if(it!=header_info.end()){
778 ss=stringstream(it->second);
779 ss>>voxelsize_x>>voxelsize_y;
780 if(ndims==3){ ss>>voxelsize_z;}
784 it=header_info.find(
"ElementSpacing");
785 _spacingX=voxelsize_x;
786 _spacingY=voxelsize_y;
787 _spacingZ=voxelsize_z;
788 if(it!=header_info.end()){
789 ss=stringstream(it->second);
790 ss>>_spacingX>>_spacingY;
791 if(ndims==3){ ss>>_spacingZ;}
800 it=header_info.find(
"Position");
801 if(it==header_info.end()){
802 it=header_info.find(
"Origin");
804 if(it==header_info.end()){
805 it=header_info.find(
"Offset");
807 if(it!=header_info.end()){
812 ss=stringstream(it->second);
813 ss>>origin_x>>origin_y;
814 if(ndims==3){ ss>>origin_z;}
816 metadata.SetDouble(
"ORIGIN_X",origin_x);
817 metadata.SetDouble(
"ORIGIN_Y",origin_y);
818 metadata.SetDouble(
"ORIGIN_Z",origin_z);
821 long nelem=dimx*dimy*dimz;
822 Init(dimx,dimy,dimz,_spacingX,_spacingY,_spacingZ);
861 string Datafile_fullpath=
"";
863 size_t found=filename.find_last_of(
"/\\");
864 if(found!=string::npos){
865 Datafile_fullpath.append(filename.substr(0,found+1));
868 Datafile_fullpath.append(datafile);
872 f.open(Datafile_fullpath,ios::binary | ios::ate);
877 int fileSize=f.tellg();
880 f.open(Datafile_fullpath,ios::in|ios::binary);
883 cout << Datafile_fullpath << endl;
887 if(typestring==
"MET_UCHAR"){
888 uint8_t* tmpdata=
new uint8_t[nelem];
889 if(BinaryDataByteOrderMSB){
896 #pragma omp parallel for
897 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
899 metadata.SetString(
"ORIGINAL_TYPE",
"UINT8");
904 else if(typestring==
"MET_CHAR"){
905 int8_t* tmpdata=
new int8_t[nelem];
906 if(BinaryDataByteOrderMSB){
913 #pragma omp parallel for
914 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
916 metadata.SetString(
"ORIGINAL_TYPE",
"INT8");
920 else if(typestring==
"MET_USHORT"){
921 uint16_t* tmpdata=
new uint16_t[nelem];
922 if(BinaryDataByteOrderMSB){
929 #pragma omp parallel for
930 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
932 metadata.SetString(
"ORIGINAL_TYPE",
"UINT16");
936 else if(typestring==
"MET_SHORT"){
937 int16_t* tmpdata=
new int16_t[nelem];
938 if(BinaryDataByteOrderMSB){
945 #pragma omp parallel for
946 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
948 metadata.SetString(
"ORIGINAL_TYPE",
"INT16");
951 else if(typestring==
"MET_UINT"){
952 uint32_t* tmpdata=
new uint32_t[nelem];
953 if(BinaryDataByteOrderMSB){
960 #pragma omp parallel for
961 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
963 metadata.SetString(
"ORIGINAL_TYPE",
"UINT32");
967 else if(typestring==
"MET_INT"){
968 int32_t* tmpdata=
new int32_t[nelem];
969 if(BinaryDataByteOrderMSB){
976 #pragma omp parallel for
977 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
979 metadata.SetString(
"ORIGINAL_TYPE",
"INT32");
983 else if(typestring==
"MET_FLOAT"){
984 float* tmpdata=
new float[nelem];
985 if(BinaryDataByteOrderMSB){
992 #pragma omp parallel for
993 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
995 metadata.SetString(
"ORIGINAL_TYPE",
"FLOAT");
998 else if(typestring==
"MET_DOUBLE"){
999 double* tmpdata=
new double[nelem];
1000 if(BinaryDataByteOrderMSB){
1007 #pragma omp parallel for
1008 for(
long i=0; i<nelem; i++){ data[i]=T(tmpdata[i]); }
1010 metadata.SetString(
"ORIGINAL_TYPE",
"DOUBLE");
1020 it=header_info.find(
"Comment");
1021 if(it!=header_info.end()){
1022 metadata.SetString(it->first,it->second);
1026 it=header_info.find(
"Name");
1027 if(it!=header_info.end()){
1028 metadata.SetString(it->first,it->second);
1032 it=header_info.find(
"ID");
1033 if(it!=header_info.end()){
1034 metadata.SetString(it->first,it->second);
1046 it=header_info.find(
"Orientation");
1047 if(it!=header_info.end()){
1048 metadata.SetString(it->first,it->second);
1052 it=header_info.find(
"AnatomicalOrientation");
1053 if(it!=header_info.end()){
1054 metadata.SetString(it->first,it->second);
1058 it=header_info.find(
"Modality");
1059 if(it!=header_info.end()){
1060 metadata.SetString(it->first,it->second);
1084 string datatype=GetDataType();
1085 string metaimage_datatype=
"";
1086 if(datatype==
"UINT8"){metaimage_datatype=
"MET_UCHAR";}
1087 else if(datatype==
"UINT16"){metaimage_datatype=
"MET_USHORT";}
1088 else if(datatype==
"UINT32"){metaimage_datatype=
"MET_UINT";}
1089 else if(datatype==
"INT8"){metaimage_datatype=
"MET_CHAR";}
1090 else if(datatype==
"INT16"){metaimage_datatype=
"MET_SHORT";}
1091 else if(datatype==
"INT32"){metaimage_datatype=
"MET_INT";}
1092 else if(datatype==
"float"){metaimage_datatype=
"MET_FLOAT";}
1093 else if(datatype==
"double"){metaimage_datatype=
"MET_DOUBLE";}
1096 outfile.open(headerfile_name,ofstream::out|ofstream::binary);
1097 if(!outfile) {
return false; }
1098 outfile<<
"ObjectType = Image"<<endl;
1101 if( GetDepth()==1){ ndims=2; }
1102 outfile<<
"NDims = "<<ndims<<endl;
1104 outfile<<
"DimSize = "<<GetWidth()<<
" "<<GetHeight();
1105 if(ndims==3){outfile<<
" "<<GetDepth();}
1108 outfile<<
"HeaderSize = 0"<<endl;
1109 outfile<<
"ElementSize = 1 1 1"<<endl;
1110 outfile<<
"ElementSpacing = "<<GetSpacingX()<<
" "<<GetSpacingY()<<
" "<<GetSpacingZ()<<
" "<<endl;
1112 double origin_x=metadata.GetDouble(
"ORIGIN_X",0);
1113 double origin_y=metadata.GetDouble(
"ORIGIN_Y",0);
1114 double origin_z=metadata.GetDouble(
"ORIGIN_Z",0);
1115 outfile<<
"Origin = "<<origin_x<<
" "<<origin_y<<
" "<<origin_z<<
" "<<endl;
1118 outfile<<
"ElementType = "<<metaimage_datatype<<endl;
1121 outfile<<
"ElementByteOrderMSB = True"<<endl;
1124 outfile<<
"ElementByteOrderMSB = False"<<endl;
1127 outfile<<
"ElementDataFile = "<<relative_datafile<<endl;
1131 long nelem= _width*_height*_depth;
1134 outfile.open(datafile,ofstream::out|ofstream::binary);
1135 outfile.write((
char*)data, nelem*
sizeof(T));
1149 f.open( filename.c_str(), std::ios::in | std::ios::binary );
1150 if( !f ) {
return false; }
1153 int32_t sizeof_hdr=0;
1154 f.read((
char*)&sizeof_hdr,4);
1155 if(sizeof_hdr!=348){
1165 f.read((
char*)&dim[0],16);
1167 if(!(dim[0]>0 && dim[0]<5)){
1174 if(dim[0]>=2){ dimy=dim[2]; }
1176 if(dim[0]>=3){ dimz=dim[3]; }
1178 if(dim[0]>=4){ nchannels=dim[4]; }
1183 f.read((
char*)&datatype,2);
1192 f.read((
char*)&pixdim[0],32);
1193 _spacingX=pixdim[1];
1194 _spacingY=pixdim[2];
1195 _spacingZ=pixdim[3];
1199 f.read((
char*)&vox_offset,4);
1202 f.seekg(std::streamoff(vox_offset));
1204 long nelem=dimx*dimy*dimz;
1205 Init(dimx,dimy,dimz);
1216 {uint8_t* tmpdata=
new uint8_t[nelem];
1217 metadata.SetString(
"ORIGINAL_TYPE",
"UINT8");
1218 f.read((
char*)tmpdata, nelem*
sizeof(uint8_t));
1219 #pragma omp parallel for
1220 for(
long i=0; i<nelem; i++){
1221 data[i]=T(tmpdata[i]);
1227 {int16_t* tmpdata=
new int16_t[nelem];
1228 metadata.SetString(
"ORIGINAL_TYPE",
"INT16");
1229 f.read((
char*)tmpdata, nelem*
sizeof(int16_t));
1230 #pragma omp parallel for
1231 for(
long i=0; i<nelem; i++){
1232 data[i]=T(tmpdata[i]);
1238 {int32_t* tmpdata=
new int32_t[nelem];
1239 metadata.SetString(
"ORIGINAL_TYPE",
"INT32");
1240 f.read((
char*)tmpdata, nelem*
sizeof(int32_t));
1241 #pragma omp parallel for
1242 for(
long i=0; i<nelem; i++){
1243 data[i]=T(tmpdata[i]);
1249 {
float* tmpdata=
new float[nelem];
1250 metadata.SetString(
"ORIGINAL_TYPE",
"FLOAT");
1251 f.read((
char*)tmpdata, nelem*
sizeof(
float));
1252 #pragma omp parallel for
1253 for(
long i=0; i<nelem; i++){
1254 data[i]=T(tmpdata[i]);
1264 {
double* tmpdata=
new double[nelem];
1265 metadata.SetString(
"ORIGINAL_TYPE",
"DOUBLE");
1266 f.read((
char*)tmpdata, nelem*
sizeof(
double));
1267 #pragma omp parallel for
1268 for(
long i=0; i<nelem; i++){
1269 data[i]=T(tmpdata[i]);
1275 {int8_t* tmpdata=
new int8_t[nelem];
1276 metadata.SetString(
"ORIGINAL_TYPE",
"INT8");
1277 f.read((
char*)tmpdata, nelem*
sizeof(int8_t));
1278 #pragma omp parallel for
1279 for(
long i=0; i<nelem; i++){
1280 data[i]=T(tmpdata[i]);
1286 {uint16_t* tmpdata=
new uint16_t[nelem];
1287 metadata.SetString(
"ORIGINAL_TYPE",
"UINT16");
1288 f.read((
char*)tmpdata, nelem*
sizeof(uint16_t));
1289 #pragma omp parallel for
1290 for(
long i=0; i<nelem; i++){
1291 data[i]=T(tmpdata[i]);
1297 {uint32_t* tmpdata=
new uint32_t[nelem];
1298 metadata.SetString(
"ORIGINAL_TYPE",
"UINT32");
1299 f.read((
char*)tmpdata, nelem*
sizeof(uint32_t));
1300 #pragma omp parallel for
1301 for(
long i=0; i<nelem; i++){
1302 data[i]=T(tmpdata[i]);
1320 outfile.open(filename);
1321 if( !outfile ) {
return false; }
1324 int32_t sizeof_hdr=348;
1325 outfile.write((
char*)&sizeof_hdr,4);
1339 outfile.write((
char*)&dim[0],16);
1344 string datatype_name=GetDataType();
1346 if(datatype_name==
"UINT8"){datatype=2;}
1347 else if(datatype_name==
"INT8"){datatype=256;}
1348 else if(datatype_name==
"UINT16"){datatype=512;}
1349 else if(datatype_name==
"INT16"){datatype=4;}
1350 else if(datatype_name==
"UINT32"){datatype=768;}
1351 else if(datatype_name==
"INT32"){datatype=8;}
1352 else if(datatype_name==
"float"){datatype=16;}
1353 else if(datatype_name==
"double"){datatype=64;}
1356 outfile.write((
char*)&datatype,2);
1362 pixdim[1]=_spacingX;
1363 pixdim[2]=_spacingY;
1364 pixdim[3]=_spacingZ;
1369 outfile.write((
char*)&pixdim[0],32);
1372 float vox_offset=352;
1374 outfile.write((
char*)&vox_offset,4);
1377 char magic[4]={
'n',
'+',
'1',
'\0'};
1379 outfile.write((
char*)&magic[0],4);
1384 outfile.write((
char*)data,
sizeof(T)*NumberOfElements());
1393 #ifdef grid3_use_gdcm
1395 gdcm::ImageReader ir;
1396 ir.SetFileName(filename.c_str());
1397 if(!ir.Read()){
return false; }
1399 const gdcm::Image &image = ir.GetImage();
1401 const unsigned int* dim = image.GetDimensions();
1402 Init(dim[0], dim[1], dim[2], image.GetSpacing(0), image.GetSpacing(1), image.GetSpacing(2));
1404 const double* origin =image.GetOrigin();
1406 metadata.SetDouble(
"ORIGIN_X",origin[0]);
1407 metadata.SetDouble(
"ORIGIN_Y",origin[1]);
1408 metadata.SetDouble(
"ORIGIN_Z",origin[2]);
1411 gdcm::PixelFormat pf= image.GetPixelFormat();
1412 if(pf.GetSamplesPerPixel()!=1){
return false;}
1414 long nelem=NumberOfElements();
1415 string type= string(pf.GetScalarTypeAsString());
1434 uint8_t* tmpdata=
new uint8_t[nelem];
1435 metadata.SetString(
"ORIGINAL_TYPE",
"UINT8");
1436 if( !image.GetBuffer( (
char*)tmpdata) ){
1440 #pragma omp parallel for
1441 for(
long i=0; i<nelem; i++){
1442 data[i]=T(tmpdata[i]);
1446 else if(type==
"UINT16"){
1447 uint16_t* tmpdata =
new uint16_t[nelem];
1448 metadata.SetString(
"ORIGINAL_TYPE",
"UINT16");
1449 if (!image.GetBuffer((
char*)tmpdata)){
1453 #pragma omp parallel for
1454 for (
long i = 0; i<nelem; i++){
1455 data[i] = T(tmpdata[i]);
1460 else if (type ==
"UINT32"){
1461 uint32_t* tmpdata =
new uint32_t[nelem];
1462 metadata.SetString(
"ORIGINAL_TYPE",
"UINT32");
1463 if (!image.GetBuffer((
char*)tmpdata)){
1467 #pragma omp parallel for
1468 for (
long i = 0; i<nelem; i++){
1469 data[i] = T(tmpdata[i]);
1473 else if (type ==
"INT8"){
1474 int8_t* tmpdata =
new int8_t[nelem];
1475 metadata.SetString(
"ORIGINAL_TYPE",
"INT8");
1476 if (!image.GetBuffer((
char*)tmpdata)){
1480 #pragma omp parallel for
1481 for (
long i = 0; i<nelem; i++){
1482 data[i] = T(tmpdata[i]);
1487 else if (type ==
"INT16"){
1488 int16_t* tmpdata =
new int16_t[nelem];
1489 metadata.SetString(
"ORIGINAL_TYPE",
"INT16");
1490 if (!image.GetBuffer((
char*)tmpdata)){
1494 #pragma omp parallel for
1495 for (
long i = 0; i<nelem; i++){
1496 data[i] = T(tmpdata[i]);
1500 else if (type ==
"INT32"){
1501 metadata.SetString(
"ORIGINAL_TYPE",
"INT32");
1502 int32_t* tmpdata =
new int32_t[nelem];
1503 if (!image.GetBuffer((
char*)tmpdata)){
1507 #pragma omp parallel for
1508 for (
long i = 0; i<nelem; i++){
1509 data[i] = T(tmpdata[i]);
1513 else if (type ==
"FLOAT16"){
1516 else if (type ==
"FLOAT32"){
1518 if (CHAR_BIT*
sizeof(
float) != 32){
1521 float* tmpdata =
new float[nelem];
1522 metadata.SetString(
"ORIGINAL_TYPE",
"FLOAT");
1523 if (!image.GetBuffer((
char*)tmpdata)){
1527 #pragma omp parallel for
1528 for (
long i = 0; i<nelem; i++){
1529 data[i] = T(tmpdata[i]);
1533 else if (type ==
"FLOAT64"){
1535 if (CHAR_BIT*
sizeof(
double) != 64){
1538 double* tmpdata =
new double[nelem];
1539 metadata.SetString(
"ORIGINAL_TYPE",
"DOUBLE");
1540 if (!image.GetBuffer((
char*)tmpdata)){
1544 #pragma omp parallel for
1545 for (
long i = 0; i<nelem; i++){
1546 data[i] = T(tmpdata[i]);
BoundaryMode
Options for reading grid data outside the grid bounds.
Definition: Grid3.h:51
std::vector< T > data
Definition: Grid3.h:59
bool WriteVTK(std::string filename)
Save grid in VTK legacy format.
Definition: Grid3_i.h:603
void ReadFromLittleEndian(char *ptr, int nelem, int bytesPerElem, std::istream &is)
Definition: Endian.h:70
T z
Definition: Vector3.h:13
int NumberOfElements()
Return the total number of grid points.
Definition: Grid3_i.h:96
Class representing 3D vectors, with elements of type T.
Definition: Vector3.h:9
T * Data()
Returns a pointer to the beginning of the internal data array. The called object owns the pointer...
Definition: Grid3_i.h:101
bool ReadVTK(std::string filename)
Construct grid from an image file in VTK legacy format.
Definition: Grid3_i.h:393
void ReadFromBigEndian(char *ptr, int nelem, int bytesPerElem, std::istream &is)
Definition: Endian.h:39
void Fill(T value)
Set all points in the grid to a specified value.
Definition: Grid3_i.h:301
std::string GetDataType()
Return a string representation of the voxel type.
Definition: Grid3_i.h:136
Grid3()
Default constructor.
Definition: Grid3.h:67
Definition: FileName.h:26
T & operator[](long index)
Access to the grid value at a specified index. The passed index is assumed to be within the grid boun...
Definition: Grid3_i.h:172
T GetMaxValue()
Return the maximum value in the grid.
Definition: Grid3_i.h:306
Vec3d GridToWorld(Vec3d p)
Transform point p from grid coordinates to world coordinates.
Definition: Grid3_i.h:87
string ToLower(string s)
Definition: Grid3_i.h:18
void Get6Neighbors(long index, vector< long > &n)
Fill the vector n with the indices of the 6-neighbors of the gridpoint at the specified index...
Definition: Grid3_i.h:241
Vec3< int > Vec3i
Definition: Vector3.h:150
T x
Definition: Vector3.h:11
void SetSpacing(double x_spacing, double y_spacing, double z_spacing)
Set spacing.
Definition: Grid3_i.h:63
void Init(int width, int height, int depth=1, double x_spacing=1.0, double y_spacing=1.0, double z_spacing=1.0)
Initialize the grid with given dimensions. Any previous data in the grid will be overwritten.
Definition: Grid3_i.h:51
void CopyDataFrom(Grid3< T > &src)
Replace the current grid data with that of another grid. If the two grids do not have the same Number...
Definition: Grid3_i.h:359
T & operator()(int x, int y, int z=0, bool Safe=false, BoundaryMode mode=REPLICATE)
Access to the grid value at a specified coordinate. If "Safe" is false, the passed coordinate is assu...
Definition: Grid3_i.h:177
Class representing a 3D image grid, each voxel being of type T.
Definition: Grid3.h:56
bool ReadMetaImage(std::string filename)
Construct grid from an image file in MetaImage format.
Definition: Grid3_i.h:677
Vec3d WorldToGrid(Vec3d p)
Transform point p from world coordinates to grid coordinates.
Definition: Grid3_i.h:77
void SwapContentsWith(Grid3< T > &v)
Swap the contents of two grids. If the two grids do not have the same NumberOfElements(), nothing happens.
Definition: Grid3_i.h:366
int Offset(int x, int y, int z)
Returns the index of a specified coordinate in the internal data array.
Definition: Grid3_i.h:142
bool WriteMetaImage(std::string filename)
Save grid in MetaImage format.
Definition: Grid3_i.h:1068
T & Background()
Access to the the background value of the grid. This is the value that is returned when trying to rea...
Definition: Grid3_i.h:167
std::string GetAsString()
Definition: FileName.h:43
bool IsBigEndian()
Definition: Endian.h:16
T LinearAt(double x, double y, double z, BoundaryMode mode=REPLICATE)
Get the grid value at a specified coordinate using trilinear interpolation. The "mode" parameter dete...
Definition: Grid3_i.h:213
bool WriteNIFTI(std::string filename)
Save grid in NIFTI format.
Definition: Grid3_i.h:1317
T GetMinValue()
Return the minimum value in the grid.
Definition: Grid3_i.h:311
void WriteToBigEndian(char *ptr, int nelem, int bytesPerElem, std::ostream &os)
Definition: Endian.h:50
bool IsInside(int x, int y, int z)
Test if the given grid point is within the bounds of the grid.
Definition: Grid3_i.h:157
void ClearPath()
Definition: FileName.h:61
bool CreateFromFile(std::string filename)
Construct grid from an image file.
Definition: Grid3_i.h:373
Vec3i Coordinate(long index)
Returns the coordinate corresponding to a specified index in the internal data array.
Definition: Grid3_i.h:147
bool ReadDICOM(std::string filename)
Construct grid from an image file in DICOM format.
Definition: Grid3_i.h:1392
T y
Definition: Vector3.h:12
void SetExtension(std::string s)
Definition: FileName.h:59
bool ReadNIFTI(std::string filename)
Construct grid from an image file in NIFTI format.
Definition: Grid3_i.h:1143