]> Creatis software - creaVtk.git/blob - bbtk_creaVtk_PKG/src/bbcreaVtkReadMHDPlane.cxx
#3061 creaVtk Bug New Immediate - ReadMHDPlane box not reading MET_INT types files
[creaVtk.git] / bbtk_creaVtk_PKG / src / bbcreaVtkReadMHDPlane.cxx
1 //===== 
2 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
3 //===== 
4 #include "bbcreaVtkReadMHDPlane.h"
5 #include "bbcreaVtkPackage.h"
6
7 #include "stdio.h"
8
9
10 #define _LARGEFILE64_SOURCE
11
12 #include <stdio.h>
13 #include <sys/types.h>
14
15 #if defined(_WIN32)
16 #else
17 #include <unistd.h>
18 #endif // defined(_WIN32)
19
20 #include <stdlib.h>
21 #include <errno.h>
22 #include <string.h>
23 #include <sys/stat.h>
24 #include <fcntl.h>
25
26 namespace bbcreaVtk
27 {
28
29 BBTK_ADD_BLACK_BOX_TO_PACKAGE(creaVtk,ReadMHDPlane)
30 BBTK_BLACK_BOX_IMPLEMENTATION(ReadMHDPlane,bbtk::AtomicBlackBox);
31 //===== 
32 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
33 //===== 
34
35
36 vtkImageData* ReadMHDPlane::CreateDefaultImage()
37 {
38         int i;
39         int sizeX, sizeY, sizeZ;
40         sizeX = 200;
41         sizeY = sizeX;
42         sizeZ = 1;
43         vtkImageData *newImage = vtkImageData::New();
44         newImage->Initialize();
45         newImage->SetScalarTypeToUnsignedChar();
46         newImage->SetSpacing( 1,1,1 );
47         newImage->SetDimensions(  sizeX,sizeY,sizeZ );
48         newImage->SetWholeExtent(0,  sizeX-1,0,sizeY-1,0,sizeZ-1 );
49         newImage->SetExtent(0,  sizeX-1,0,sizeY-1,0,sizeZ-1 );
50         newImage->SetNumberOfScalarComponents(1);
51         newImage->AllocateScalars();
52         newImage->Update();
53         memset ( (void*)newImage->GetScalarPointer(), 0, sizeX*sizeY*1 );
54         for (i=0; i<sizeX; i++)
55         {
56                 newImage->SetScalarComponentFromDouble(i,i,0, 0, 255 );
57                 newImage->SetScalarComponentFromDouble(i,sizeY-1-i,0, 0, 255 );
58         } // for i
59         return newImage;
60
61
62
63 void ReadMHDPlane::ReadNormalMHD()
64 {
65 // THE MAIN PROCESSING METHOD BODY
66 //   Here we simply set the input 'In' value to the output 'Out'
67 //   And print out the output value
68 // INPUT/OUTPUT ACCESSORS ARE OF THE FORM :
69 //    void bbSet{Input|Output}NAME(const TYPE&)
70 //    const TYPE& bbGet{Input|Output}NAME() const 
71 //    Where :
72 //    * NAME is the name of the input/output
73 //      (the one provided in the attribute 'name' of the tag 'input')
74 //    * TYPE is the C++ type of the input/output
75 //      (the one provided in the attribute 'type' of the tag 'input')
76
77 //    bbSetOutputOut( bbGetInputIn() );
78 //    std::cout << "Output value = " <<bbGetOutputOut() << std::endl;
79
80         std::string inputfilename;
81         int slice;
82         int width;
83         long long dataSize;
84
85
86         if (bbGetInputDirectionPlane()=="XY")
87         {
88                 inputfilename=bbGetInputFileName();
89         } // if XY
90
91
92         if ((bbGetInputDirectionPlane()=="YZ") || (bbGetInputDirectionPlane()=="ZX"))
93         {
94                 //-- Split FileName
95                 std::size_t found               =       bbGetInputFileName().find_last_of("/\\");
96                 std::string path                =       bbGetInputFileName().substr(0,found+1);
97                 std::string filename    =       bbGetInputFileName().substr(found+1);
98 #ifdef _WIN32
99                 path=path+"YZ_ZX\\";
100 #else
101                 path=path+"YZ_ZX/";
102 #endif
103                 if (bbGetInputDirectionPlane()=="YZ") 
104                 { 
105                         inputfilename = path+filename+"_YZ.mhd";  
106                 } // if YZ
107
108                 if (bbGetInputDirectionPlane()=="ZX")
109                 {
110                         inputfilename = path+filename+"_ZX.mhd";
111                 } // XZ
112         } // if YZ || XZ
113
114         width = bbGetInputWidth();
115         if (width<=0 ) { width=1; }
116
117         slice = bbGetInputSlice();
118         if (slice<0 ) { slice=0; }
119
120         vtkImageData *newImage;
121         long long newHeaderSize;
122     std::string newFileName=inputfilename+"-OneSlice";
123
124     int sx,sy,sz;
125     char mystring[250];
126     char strTmp[20];
127         bool ok=true;
128     FILE *ffIn  = fopen(inputfilename.c_str(),"r");
129         if (ffIn!=NULL)
130         {
131             FILE *ffOut = fopen(newFileName.c_str(),"w");
132
133         while(!feof(ffIn))
134         {
135             strcpy(mystring,"\n");
136                         fgets(mystring,250,ffIn);
137                 if (strncmp("NDims",mystring,5)==0) 
138                         {
139                                 if (width==1) 
140                                 { 
141                                         strcpy(mystring,"NDims = 2\n"); 
142                                 } else {
143                                         strcpy(mystring,"NDims = 3\n"); 
144                                 }
145                         }
146                 if (strncmp("DimSize",mystring,6)==0) 
147                         {
148                                 sscanf(mystring,"%s %s %d %d %d",strTmp, strTmp, &sx, &sy,&sz);
149
150                                 if (width==1) 
151                                 { 
152                                         sprintf(mystring,"DimSize = %d %d\n",sx,sy);
153                                 } else {
154                                         sprintf(mystring,"DimSize = %d %d %d\n",sx,sy,width);
155                                 }
156
157                         newHeaderSize = (long long int)sx*(long long int)sy*(long long int)slice;
158                                 if (bbGetInputSlice()>=sz) {ok=false;}
159                         } // if
160                 if (strncmp("ElementType = MET_CHAR",mystring,22)==0)                   { dataSize=sizeof(char);                        }
161                 if (strncmp("ElementType = VTK_CHAR",mystring,22)==0)                   { dataSize=sizeof(char);                        }
162                 if (strncmp("ElementType = MET_UCHAR",mystring,23)==0)                  { dataSize=sizeof(unsigned char);       }
163                 if (strncmp("ElementType = VTK_UNSIGNED_CHAR",mystring,31)==0)  { dataSize=sizeof(unsigned char);       }
164                 if (strncmp("ElementType = MET_USHORT",mystring,24)==0)                 { dataSize=sizeof(unsigned short);      }
165                 if (strncmp("ElementType = VTK_UNSIGNED_SHORT",mystring,32)==0) { dataSize=sizeof(unsigned short);      }
166                 if (strncmp("ElementType = MET_SHORT",mystring,23)==0)                  { dataSize=sizeof(short);                       }
167                 if (strncmp("ElementType = VTK_SHORT",mystring,23)==0)                  { dataSize=sizeof(short);                       }
168                 if (strncmp("ElementType = MET_UINT",mystring,22)==0)                   { dataSize=sizeof(unsigned int);        }
169                 if (strncmp("ElementType = VTK_UNSIGNED_INT",mystring,30)==0)   { dataSize=sizeof(unsigned int);        }
170                 if (strncmp("ElementType = MET_INT",mystring,21)==0)                    { dataSize=sizeof(int);                         }
171                 if (strncmp("ElementType = VTK_INT",mystring,21)==0)                    { dataSize=sizeof(int);                         }
172                 if (strncmp("ElementType = MET_FLOAT",mystring,23)==0)                  { dataSize=sizeof(float);                       }
173                 if (strncmp("ElementType = VTK_FLOAT",mystring,23)==0)                  { dataSize=sizeof(float);                       }
174                 if (strncmp("ElementType = MET_LONG",mystring,22)==0)                   { dataSize=sizeof(long);                        }
175                 if (strncmp("ElementType = VTK_LONG",mystring,22)==0)                   { dataSize=sizeof(long);                        }
176                 if (strncmp("ElementType = MET_DOUBLE",mystring,24)==0)                 { dataSize=sizeof(double);                      }
177                 if (strncmp("ElementType = VTK_DOUBLE",mystring,24)==0)                 { dataSize=sizeof(double);                      }
178                         newHeaderSize=newHeaderSize*dataSize;
179
180                 if (strncmp("Offset",mystring,6)==0) {strcpy(mystring,"Offset = 0 0 0\n");}
181                 if (strncmp("HeaderSize",mystring,10)==0) {strcpy(mystring,"");}
182                 if (strncmp("ElementDataFile",mystring,15)==0) 
183                         {
184                                 fprintf(ffOut,"HeaderSize = %lld\n\n", newHeaderSize );
185                         } // if
186                 fprintf(ffOut,mystring); 
187         } // while
188                 fclose(ffIn);
189                 fclose(ffOut);
190
191                 if (ok==true)
192                 {
193                         vtkMetaImageReader * reader = vtkMetaImageReader::New();
194                         reader->SetFileName( newFileName.c_str() );
195                         reader->Update();
196                         newImage = reader->GetOutput();
197                 } // if ok
198         } else {
199                 ok=false;
200                 printf("EED ERROR: Problem openin:%s\n", inputfilename.c_str() );
201         }
202
203         if (ok==false)
204         {
205                 newImage=CreateDefaultImage();
206         } // if ok
207
208         bbSetOutputOut( newImage );
209 }
210
211
212
213
214
215 void ReadMHDPlane::Read64lseek()
216 {
217         int             slice;
218         int             width;
219
220         width = bbGetInputWidth();
221         if (width<=0 ) { width=1; }
222
223         slice = bbGetInputSlice();
224         if (slice<0 ) { slice=0; }
225
226         int             dimX=-1,dimY=-1,dimZ=-1;
227         int             dim=-1;
228         std::string formattype("VOID");
229         std::string elementdatafile("VOID");
230         float           spcX=-1,spcY=-1,spcZ=-1;
231         float           ox=-1,oy=-1,oz=-1;
232         long int        headersize=0;
233
234         vtkImageData *newImage=NULL;
235
236         char mystring[250];
237         char strTmp[30];
238         char strTmp2[30];
239         FILE *ffIn      = fopen(bbGetInputFileName().c_str(),"r");
240         long long dataSize;
241
242         if (ffIn!=NULL)
243         {
244                 newImage = vtkImageData::New();
245
246                 while(!feof(ffIn))
247                 {
248             strcpy(mystring,"\n");
249                         fgets(mystring,250,ffIn);
250                 if (strncmp("NDims",mystring,5)==0)                     { sscanf(mystring,"%s %s %d"            ,strTmp, strTmp, &dim);                                 }
251                 if (strncmp("DimSize",mystring,6)==0)                   { sscanf(mystring,"%s %s %d %d %d"      ,strTmp, strTmp, &dimX, &dimY,&dimZ);   } 
252                         if (strncmp("ElementType",mystring,11)==0)              { sscanf(mystring,"%s %s %s"            ,strTmp, strTmp, strTmp2); formattype=strTmp2;                  }
253                         if (strncmp("ElementSpacing",mystring,14)==0)   { sscanf(mystring,"%s %s %f %f %f"      ,strTmp, strTmp, &spcX,&spcY,&spcZ);    }
254                         if (strncmp("ElementSize",mystring,11)==0)              { sscanf(mystring,"%s %s %f %f %f"      ,strTmp, strTmp, &spcX,&spcY,&spcZ);    }
255                 if (strncmp("Offset",mystring,6)==0)                    { sscanf(mystring,"%s %s %f %f %f"      ,strTmp, strTmp, &ox, &oy, &oz);                }
256                 if (strncmp("HeaderSize",mystring,10)==0)               { sscanf(mystring,"%s %s %ld"           ,strTmp, strTmp, &headersize);                  }
257                         if (strncmp("ElementDataFile",mystring,15)==0)  { sscanf(mystring,"%s %s %s"            ,strTmp, strTmp, strTmp2); elementdatafile=strTmp2;             }
258
259                 if (strncmp("ElementType = MET_CHAR",mystring,22)==0)                   { newImage->SetScalarTypeToChar();                      dataSize=sizeof(char);                  }
260                 if (strncmp("ElementType = VTK_CHAR",mystring,22)==0)                   { newImage->SetScalarTypeToChar();                      dataSize=sizeof(char);                  }
261                 if (strncmp("ElementType = MET_UCHAR",mystring,23)==0)                  { newImage->SetScalarTypeToUnsignedChar();      dataSize=sizeof(unsigned char); }
262                 if (strncmp("ElementType = VTK_UNSIGNED_CHAR",mystring,31)==0)  { newImage->SetScalarTypeToUnsignedChar();      dataSize=sizeof(unsigned char); }
263                 if (strncmp("ElementType = MET_USHORT",mystring,24)==0)                 { newImage->SetScalarTypeToUnsignedShort();     dataSize=sizeof(unsigned short);}
264                 if (strncmp("ElementType = VTK_UNSIGNED_SHORT",mystring,32)==0) { newImage->SetScalarTypeToUnsignedShort();     dataSize=sizeof(unsigned short);}
265                 if (strncmp("ElementType = MET_SHORT",mystring,23)==0)                  { newImage->SetScalarTypeToShort();                     dataSize=sizeof(short);                 }
266                 if (strncmp("ElementType = VTK_SHORT",mystring,23)==0)                  { newImage->SetScalarTypeToShort();                     dataSize=sizeof(short);                 }
267                 if (strncmp("ElementType = MET_UINT",mystring,22)==0)                   { newImage->SetScalarTypeToUnsignedInt();       dataSize=sizeof(unsigned int);  }
268                 if (strncmp("ElementType = VTK_UNSIGNED_INT",mystring,30)==0)   { newImage->SetScalarTypeToUnsignedInt();       dataSize=sizeof(unsigned int);  }
269                 if (strncmp("ElementType = MET_INT",mystring,21)==0)                    { newImage->SetScalarTypeToInt();                       dataSize=sizeof(int);                   }
270                 if (strncmp("ElementType = VTK_INT",mystring,21)==0)                    { newImage->SetScalarTypeToInt();                       dataSize=sizeof(int);                   }
271                 if (strncmp("ElementType = MET_FLOAT",mystring,23)==0)                  { newImage->SetScalarTypeToFloat();                     dataSize=sizeof(float);                 }
272                 if (strncmp("ElementType = VTK_FLOAT",mystring,23)==0)                  { newImage->SetScalarTypeToFloat();                     dataSize=sizeof(float);                 }
273                 if (strncmp("ElementType = MET_LONG",mystring,22)==0)                   { newImage->SetScalarTypeToLong();                      dataSize=sizeof(long);                  }
274                 if (strncmp("ElementType = VTK_LONG",mystring,22)==0)                   { newImage->SetScalarTypeToLong();                      dataSize=sizeof(long);                  }
275                 if (strncmp("ElementType = MET_DOUBLE",mystring,24)==0)                 { newImage->SetScalarTypeToDouble();            dataSize=sizeof(double);                }
276                 if (strncmp("ElementType = VTK_DOUBLE",mystring,24)==0)                 { newImage->SetScalarTypeToDouble();            dataSize=sizeof(double);                }
277
278         } // while
279                 fclose(ffIn);
280
281                 newImage->Initialize();
282
283
284                 int             fd;
285                 long long       ret;
286                 std::size_t found;
287                 std::string filename;
288                 found                                   = bbGetInputFileName().find_last_of("/\\");
289                 filename                                = bbGetInputFileName().substr(0,found+1) + elementdatafile ;
290                 long long       pos;
291
292                 long long lsize                 = dimX*dimY*width *dataSize;
293
294
295 #if defined(_WIN32)
296                    _sopen_s( &fd, filename.c_str(), _O_RDONLY, _SH_DENYNO, 0 );
297 #else
298                 fd = open ( filename.c_str() ,  O_RDONLY|O_LARGEFILE );
299 #endif // defined(_WIN32)
300
301
302                 if (fd < 0)
303                 {
304                         printf("EED ReadMHDPlane::Read64lseek   WARNNING! raw file not exist\n");
305                         fprintf(stderr, "%s\n", strerror(errno));
306                         newImage=CreateDefaultImage();
307 //                      exit(1);
308                 }       
309                 
310
311                 if ((bbGetInputDirectionPlane()=="XY") && (fd>=0)) 
312                 { 
313                         newImage->SetSpacing( spcX,spcY,spcZ );
314                         newImage->SetDimensions(  dimX,dimY,width );
315                         newImage->SetWholeExtent(0,  dimX-1,0,dimY-1,0,width-1 );
316                         newImage->SetExtent(0,  dimX-1,0,dimY-1,0,width-1 );
317                         newImage->SetNumberOfScalarComponents(1);
318                         newImage->AllocateScalars();
319                         newImage->Update();
320                         pos = dimX*dimY*(long long)slice*dataSize;
321 #if defined(_WIN32)
322                         if (_lseeki64( fd, pos, SEEK_SET ) < 0)
323 #else
324                         if (lseek64(fd, pos, SEEK_SET) < 0) 
325 #endif // defined(_WIN32)
326                         {
327                                 printf("EED ReadMHDPlane::Read64lseek \n");
328                                 fprintf(stderr, "Failed seeking to %lld, %s\n", pos, strerror(errno));
329                                 exit(1);
330                         }
331
332                         if ((ret = read(fd, newImage->GetScalarPointer() , lsize)) < 0) 
333                         {
334                                 fprintf(stderr, "Failed reading: %s\n", strerror(errno));
335                                 exit(1);
336                         }
337                 }  // if PLANE XY
338
339                 if ((bbGetInputDirectionPlane()=="XZ") && (fd>=0))
340                 { 
341                         newImage->SetSpacing( spcX,spcZ,spcY );
342                         newImage->SetDimensions(  dimX,dimZ,width );
343                         newImage->SetWholeExtent(0,  dimX-1,0,dimZ-1,0,width-1 );
344                         newImage->SetExtent(0,  dimX-1,0,dimZ-1,0,width-1 );
345                         newImage->SetNumberOfScalarComponents(1);
346                         newImage->AllocateScalars();
347                         newImage->Update();
348                         int iWidth;
349                         for (iWidth=0;iWidth<width;iWidth++)
350                         {
351                                 copy_XZ_plane(fd,newImage,slice+iWidth,iWidth,dimX,dimY,dimZ,dataSize);
352                         }
353                 }  // if PLANE XZ
354
355                 if ((bbGetInputDirectionPlane()=="YZ") && (fd>=0))
356                 { 
357                         newImage->SetSpacing( spcY,spcZ,spcX );
358                         newImage->SetDimensions(  dimY,dimZ,width );
359                         newImage->SetWholeExtent(0,  dimY-1,0,dimZ-1,0,width-1 );
360                         newImage->SetExtent(0,  dimY-1,0,dimZ-1,0,width-1 );
361                         newImage->SetNumberOfScalarComponents(1);
362                         newImage->AllocateScalars();
363                         newImage->Update();
364                         int iWidth;
365                         for (iWidth=0;iWidth<width;iWidth++)
366                         {
367                                 copy_YZ_plane(fd,newImage,slice+iWidth,iWidth,dimX,dimY,dimZ,dataSize);
368                         }
369
370                 }  // if PLANE YZ       
371 #if defined(_WIN32)
372                 _close (fd);
373 #else
374                 close (fd);
375 #endif // defined(_WIN32)
376         } else {
377                 newImage=CreateDefaultImage();
378         } // if  ffIn
379
380         bbSetOutputOut( newImage );
381 }
382
383
384 void ReadMHDPlane::copy_YZ_plane(int fd,vtkImageData *newImage,int slice,int iWidth,int dimX,int dimY,int dimZ,int dataSize)
385 {
386         long long int j;
387         long long int i;
388         long long       ret;
389         long long       pos;
390         long long sizeBytesPlane = dimX*dimY*dataSize;
391         char *pImage;
392         for (j=0;j<dimZ;j++)
393         {
394                 for (i=0;i<dimY;i++)
395                 {
396                         pos = ((long long int)slice+iWidth + i*dimX + j*dimX*dimY)*dataSize ;
397 #if defined(_WIN32)
398                         if (_lseeki64( fd, pos , SEEK_SET ) < 0)
399 #else
400                         if (lseek64(fd, pos , SEEK_SET) < 0) 
401 #endif // defined(_WIN32)
402                         {
403                                 printf("EED ReadMHDPlane::Read64lseek \n");
404                                 fprintf(stderr, "Failed seeking to %lld, %s\n", pos, strerror(errno));
405                                 exit(1);
406                         }
407
408                         pImage=(char*)(newImage->GetScalarPointer(i, dimZ-1-j,iWidth ));
409                         if ((ret = read(fd, pImage , dataSize)) < 0) 
410                         {
411                                 fprintf(stderr, "Failed reading: %s\n", strerror(errno));
412                                 exit(1);
413                         }
414                         } // for i
415                 } // for j
416 }
417
418
419 void ReadMHDPlane::copy_XZ_plane(int fd,vtkImageData *newImage,int slice,int iWidth,int dimX,int dimY,int dimZ,int dataSize)
420 {
421         long long int j;
422         long long       ret;
423         long long       pos;
424         pos = dimX*(long long int)slice*dataSize;
425         long long sizeBytesPlane = dimX*dimY*dataSize;
426         char *pImage = (char*)( newImage->GetScalarPointer(0,0,iWidth) );
427         for (j=dimZ-1;j>=0;j--)
428         {
429 #if defined(_WIN32)
430                 if (_lseeki64( fd, pos + j*sizeBytesPlane , SEEK_SET ) < 0)
431 #else
432                 if (lseek64(fd, pos + j*sizeBytesPlane , SEEK_SET) < 0) 
433 #endif // defined(_WIN32)
434                 {
435                         printf("EED ReadMHDPlane::Read64lseek \n");
436                         fprintf(stderr, "Failed seeking to %lld, %s\n", pos, strerror(errno));
437                         exit(1);
438                 }
439
440                 if ((ret = read(fd, pImage , dimX*dataSize)) < 0) 
441                 {
442                         fprintf(stderr, "Failed reading: %s\n", strerror(errno));
443                         exit(1);
444                 }
445                 pImage = pImage+dimX*dataSize;
446         } // for j
447 }
448
449
450
451 void ReadMHDPlane::Process()
452 {
453         if (bbGetInputType()==0)
454         {
455                 ReadNormalMHD();
456         }       
457
458         if (bbGetInputType()==1)
459         {
460                 Read64lseek();
461         }       
462
463 }
464 //===== 
465 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
466 //===== 
467 void ReadMHDPlane::bbUserSetDefaultValues()
468 {
469
470 //  SET HERE THE DEFAULT INPUT/OUTPUT VALUES OF YOUR BOX 
471 //    Here we initialize the input 'In' to 0
472    bbSetInputFileName("");
473    bbSetInputSlice(0);
474    bbSetInputType(1);
475    bbSetInputWidth(1);
476    bbSetInputDirectionPlane("XY");
477   
478 }
479 //===== 
480 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
481 //===== 
482 void ReadMHDPlane::bbUserInitializeProcessing()
483 {
484
485 //  THE INITIALIZATION METHOD BODY :
486 //    Here does nothing 
487 //    but this is where you should allocate the internal/output pointers 
488 //    if any 
489
490   
491 }
492 //===== 
493 // Before editing this file, make sure it's a file of your own (i.e.: it wasn't generated from xml description; if so : your modifications will be lost)
494 //===== 
495 void ReadMHDPlane::bbUserFinalizeProcessing()
496 {
497
498 //  THE FINALIZATION METHOD BODY :
499 //    Here does nothing 
500 //    but this is where you should desallocate the internal/output pointers 
501 //    if any
502   
503 }
504 }
505 // EO namespace bbcreaVtk
506
507