A.I.M.S


spmR_d.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 /*
35  * SPM reader template functions
36  */
37 
38 #include <cstdlib>
39 #include <aims/io/spmR.h>
40 #include <aims/io/datatypecode.h>
41 #include <aims/io/defaultItemR.h>
42 #include <aims/io/spmheader.h>
43 #include <aims/utility/flip.h>
44 #include <aims/resampling/motion.h>
48 
49 namespace aims
50 {
51 
52  template< class T >
53  std::string SpmReader< T >::removeExtension(const std::string& name)
54  {
55  std::string res = name;
56  std::string ext="";
57  if ( res.length() > 4 )
58  ext = res.substr( res.length() - 4, 4 );
59  if (ext == ".hdr" || ext == ".img")
60  res = res.substr( 0, res.length() - 4 );
61  return res;
62  }
63 
64 
65  template <typename T>
67  const carto::AllocatorContext & context,
68  carto::Object options )
69  {
70  std::string name = removeExtension( _name );
71 
72  aims::SpmHeader *hdr = new aims::SpmHeader( name + ".hdr" );
73 
74  try
75  {
76  hdr->read();
77  }
78  catch( std::exception & e )
79  {
80  delete hdr;
81  throw e;
82  }
83 
84  std::string type;
85  if ( !hdr->getProperty( "disk_data_type", type ) )
87  hdr->setProperty( "data_type", carto::DataTypeCode<T>().dataType() );
88  hdr->setProperty( "bits_allocated", (int) sizeof( T ) );
89 
90  std::vector<int> dims;
91  std::vector<float> vs;
92  if( ! hdr->getProperty( "volume_dimension", dims ) )
93  throw std::logic_error( "Internal error: getProperty volume_dimension failed" );
94  if( ! hdr->getProperty( "voxel_size", vs ) )
95  throw std::logic_error( "Internal error: getProperty voxel_size failed" );
96 
97  int frame = -1, border = 0;
98  options->getProperty( "frame", frame );
99  options->getProperty( "border", border );
100 
101  int dimt = 1, tmin = 0;
102  if ( frame >= 0 )
103  {
104  if ( frame >= dims[3] )
105  {
106  throw std::domain_error( "Frame number exceeds "
107  "volume time dimension" );
108  }
109  tmin = frame;
110  }
111  else
112  dimt = dims[3];
113 
114  int bswapi = false;
115  hdr->getProperty( "byte_swapping", bswapi );
116  bool bswap = bswapi;
117  std::string filename = name + ".img";
118 
119  bool series = false;
120  std::vector<carto::Object> fnames;
121  std::string bname;
122  if ( tmin > 0 || dimt > 1 )
123  {
124  series = hdr->getProperty( "series_filenames", fnames );
125  if ( series )
126  {
127  if ( fnames.size() < (unsigned) ( dimt + tmin ) )
128  throw std::domain_error( "SPM reader: not enough filenames for "
129  "all frames" );
130  bname = carto::FileUtil::dirname( name )
132  hdr->removeProperty( "series_filenames" );
133  }
134  }
135 
136  carto::AllocatorContext al
137  ( context.accessMode(),
139  ( new carto::FileDataSource( filename, 0, carto::DataSource::Read ) ),
140  false, context.useFactor() );
141 
142  //std::cout << "alloc type: " << al.allocatorType() << std::endl;
143  AimsData<T> data( dims[0], dims[1], dims[2], dimt, border, al );
144  data.setSizeXYZT( vs[0], vs[1] , vs[2], vs[3] );
145 
146  std::ifstream is;
147  if ( !series )
148  {
149  is.open( filename.c_str(), std::ios::in | std::ios::binary );
150  if ( !is )
151  {
152  carto::io_error::launchErrnoExcept( filename );
153  }
154  is.unsetf( std::ios::skipws );
155  if ( type == "S8")
156  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
157  * data.dimX() * sizeof( int8_t ), std::ios::cur );
158  else if ( type == "U8")
159  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
160  * data.dimX() * sizeof( uint8_t ), std::ios::cur );
161  else if ( type =="S16" )
162  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
163  * data.dimX() * sizeof( int16_t ), std::ios::cur );
164  else if ( type =="U16" )
165  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
166  * data.dimX() * sizeof( uint16_t ), std::ios::cur );
167  else if ( type =="S32" )
168  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
169  * data.dimX() * sizeof( int32_t ), std::ios::cur );
170  else if ( type =="U32" )
171  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
172  * data.dimX() * sizeof( uint32_t ), std::ios::cur );
173  else if ( type == "FLOAT" )
174  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
175  * data.dimX() * sizeof( float ), std::ios::cur );
176  else if ( type == "DOUBLE" )
177  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
178  * data.dimX() * sizeof( double ), std::ios::cur );
179  else
180  is.seekg( (std::ifstream::pos_type) tmin * data.dimZ() * data.dimY()
181  * data.dimX() * sizeof( T ), std::ios::cur );
182  }
183 
184  std::vector<T> buffer ( dims[0] );
185  std::string fnm;
186 
187  std::vector<float> vstom;
188  Motion stom;
189  if( hdr->getProperty( "storage_to_memory", vstom ) )
190  stom = Motion( vstom );
191  else
192  {
193  // standard orientation
194  stom.rotation()(1,1) = -1;
195  stom.rotation()(2,2) = -1;
196  stom.translation()[1] = data.dimY() - 1;
197  stom.translation()[2] = data.dimZ() - 1;
198  bool radio = false;
199  try
200  {
201  radio = (bool) hdr->getProperty( "spm_radio_convention" )->getScalar();
202  }
203  catch( std::exception & )
204  {
205  }
206  if( !radio )
207  {
208  stom.rotation()(0,0) = -1;
209  stom.translation()[0] = data.dimX() - 1;
210  }
211  }
212 
213  for ( int t=0; t<dimt; ++t )
214  {
215  if ( series )
216  {
217  fnm = bname
218  + carto::FileUtil::basename( fnames[t+tmin]->getString() );
219  is.open( ( fnm ).c_str(), std::ios::in | std::ios::binary );
220  if ( !is )
221  {
222  carto::io_error::launchErrnoExcept( fnm );
223  }
224  is.unsetf( std::ios::skipws );
225  }
226  hdr->setProperty( "scale_factor_applied", false ) ;
227  readFrame( is, data, t, type, hdr, bswap, stom );
228  if ( series )
229  is.close();
230  }
231 
232  if ( !series )
233  is.close();
234 
235  thing = data;
236  thing.setHeader( hdr );
237  }
238 
239 
240  template<typename T> void
241  SpmReader<T>::readFrame( std::ifstream & is, AimsData<T> & data, int t,
242  const std::string &, SpmHeader*, bool bswap,
243  const Motion & stom )
244  {
245  DefaultItemReader<T> itemR ;
246  ItemReader<T> *ir
247  = itemR.reader( "binar", bswap );
248 
249  Point3df incf = stom.transform( Point3df( 1, 0, 0 ) )
250  - stom.transform( Point3df( 0, 0, 0 ) );
251  Point4d inc = Point4d( int16_t( rint( incf[0] ) ),
252  int16_t( rint( incf[1] ) ),
253  int16_t( rint( incf[2] ) ), 0 );
254  Point3df d0f;
255  Point4d d0;
256  Motion m2s = stom.inverse();
257  Point3df df = m2s.transform( Point3df( data.dimX(), data.dimY(),
258  data.dimZ() ) )
259  - m2s.transform( Point3df( 0, 0, 0 ) );
260  Point4d tdims = Point4d( short( rint( fabs( df[0] ) ) ),
261  short( rint( fabs( df[1] ) ) ),
262  short( rint( fabs( df[2] ) ) ), (short) data.dimT() );
263  int x, dx = tdims[0];
264  std::vector<T> buffer( dx );
265 
266  for ( int z=0;z<tdims[2];++z)
267  for ( int y=0;y<tdims[1];++y)
268  {
269  d0f = stom.transform( Point3df( 0, y, z ) );
270  d0 = Point4d( int16_t( rint( d0f[0] ) ), int16_t( rint( d0f[1] ) ),
271  int16_t( rint( d0f[2] ) ), t );
272  if( inc == Point4d( 1, 0, 0, 0 ) )
273  ir->read( is, &data( d0 ), dx );
274  else
275  {
276  ir->read( is, &buffer[0], dx );
277  for ( x=0; x<dx; ++x, d0+=inc )
278  data( d0 ) = buffer[ x ];
279  }
280  }
281  delete ir ;
282  }
283 
284 
285  template <typename T> template<typename U>
286  void
287  SpmReader<T>::readScaledFrame( std::ifstream & is, AimsData<T>& data,
288  int t, SpmHeader *hdr, bool bswap,
289  const Motion & stom )
290  {
291  float scaleFactor = 1;
292  bool scaleFactorApplied = false;
293 
294  hdr->getProperty( "scale_factor", scaleFactor );
295  hdr->getProperty( "scale_factor_applied", scaleFactorApplied );
296  if ( scaleFactorApplied || scaleFactor == 0. )
297  scaleFactor = 1. ;
298 
299  DefaultItemReader<U> itemR ;
300  ItemReader<U> *ir
301  = itemR.reader( "binar", bswap );
302 
303  Point3df incf = stom.transform( Point3df( 1, 0, 0 ) )
304  - stom.transform( Point3df( 0, 0, 0 ) );
305  Point4d inc = Point4d( int16_t( rint( incf[0] ) ),
306  int16_t( rint( incf[1] ) ),
307  int16_t( rint( incf[2] ) ), 0 );
308  Point3df d0f;
309  Point4d d0;
310  Motion m2s = stom.inverse();
311  Point3df df = m2s.transform( Point3df( data.dimX(), data.dimY(),
312  data.dimZ() ) )
313  - m2s.transform( Point3df( 0, 0, 0 ) );
314  Point4d tdims = Point4d( short( rint( fabs( df[0] ) ) ),
315  short( rint( fabs( df[1] ) ) ),
316  short( rint( fabs( df[2] ) ) ), data.dimT() );
317  int x, dx = tdims[0];
318  std::vector<U> buffer( dx );
319 
320  for ( int z=0;z<tdims[2];++z)
321  for ( int y=0;y<tdims[1];++y)
322  {
323  d0f = stom.transform( Point3df( 0, y, z ) );
324  d0 = Point4d( int16_t( rint( d0f[0] ) ), int16_t( rint( d0f[1] ) ),
325  int16_t( rint( d0f[2] ) ), t );
326  ir->read( is, &buffer[0] , dx );
327  for ( x=0; x<dx; ++x, d0+=inc )
328  data( d0 ) = (T) ( scaleFactor * buffer[ x ] );
329  }
330  delete ir;
331  hdr->setProperty( "scale_factor_applied", true ) ;
332  }
333 
334 
335  template<> void
336  SpmReader<float>::readFrame( std::ifstream & is, AimsData<float> & data,
337  int t, const std::string & type,
338  SpmHeader* hdr, bool bswap,
339  const Motion & stom )
340  {
341  if ( type == "S8")
342  readScaledFrame<int8_t>( is, data, t, hdr, bswap, stom );
343  else if ( type == "U8")
344  readScaledFrame<uint8_t>( is, data, t, hdr, bswap, stom );
345  else if ( type =="S16" )
346  readScaledFrame<int16_t>( is, data, t, hdr, bswap, stom );
347  else if ( type =="U16" )
348  readScaledFrame<uint16_t>( is, data, t, hdr, bswap, stom );
349  else if ( type =="S32" )
350  readScaledFrame<int32_t>( is, data, t, hdr, bswap, stom );
351  else if ( type =="U32" )
352  readScaledFrame<uint32_t>( is, data, t, hdr, bswap, stom );
353  else if ( type =="FLOAT" )
354  readScaledFrame<float>( is, data, t, hdr, bswap, stom );
355  else if ( type =="DOUBLE" )
356  readScaledFrame<double>( is, data, t, hdr, bswap, stom );
357  else
358  std::cerr << "SPM read not implemented: " << type << " as FLOAT"
359  << std::endl;
360  }
361 
362  template<> void
363  SpmReader<double>::readFrame( std::ifstream & is, AimsData<double>& data,
364  int t, const std::string & type,
365  SpmHeader* hdr, bool bswap,
366  const Motion & stom )
367  {
368  if ( type == "S8")
369  readScaledFrame<int8_t>( is, data, t, hdr, bswap, stom );
370  else if ( type == "U8")
371  readScaledFrame<uint8_t>( is, data, t, hdr, bswap, stom );
372  else if ( type =="S16" )
373  readScaledFrame<int16_t>( is, data, t, hdr, bswap, stom );
374  else if ( type =="U16" )
375  readScaledFrame<uint16_t>( is, data, t, hdr, bswap, stom );
376  else if ( type =="S32" )
377  readScaledFrame<int32_t>( is, data, t, hdr, bswap, stom );
378  else if ( type =="U32" )
379  readScaledFrame<uint32_t>( is, data, t, hdr, bswap, stom );
380  else if ( type =="FLOAT" )
381  readScaledFrame<float>( is, data, t, hdr, bswap, stom );
382  else if ( type =="DOUBLE" )
383  readScaledFrame<double>( is, data, t, hdr, bswap, stom );
384  else
385  std::cerr << "SPM read not implemented: " << type << " as DOUBLE"
386  << std::endl;
387  }
388 
389 }
390 
SPM Header class.
Definition: spmheader.h:53
std::string dataType()
virtual bool getProperty(const std::string &key, Object &value) const =0
int dimX() const
static std::string basename(const std::string &)
int dimT() const
Default low-levels readers.
Definition: defaultItemR.h:55
virtual bool removeProperty(const std::string &)
void setHeader(aims::Header *hdr)
static char separator()
virtual ItemReader< T > * reader(const std::string &openmode="binar", bool bswap=false) const
Definition: defaultItemR.h:154
AimsVector< int16_t, 4 > Point4d
Definition: vector.h:213
AimsVector< float, 3 > Point3df
Definition: vector.h:237
Point3dd transform(double x, double y, double z) const
aims::AffineTransformation3d Motion
Definition: graphmanip.h:48
static std::string dirname(const std::string &)
void read(AimsData< T > &thing, const carto::AllocatorContext &context, carto::Object options)
Read the data with "name" file name from disk.
Definition: spmR_d.h:66
virtual ItemReader< T > * reader(const std::string &openmode="binar", bool bswap=false) const =0
virtual void read(std::istream &is, T &item) const
Definition: itemR.h:103
The template class for SPM read operation.
Definition: spmR.h:57
virtual bool getProperty(const std::string &, Object &) const
int dimZ() const
virtual void setProperty(const std::string &, Object)
int dimY() const
std::string removeExtension(const std::string &name)
Return a name without .hdr or .img extension.
Definition: spmR_d.h:53
AffineTransformation3d inverse() const
Low-level "small item" reader, used by higher-level file readers.
Definition: itemR.h:98
Affine 3D transformation.