A.I.M.S


spmW_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 #ifndef AIMS_IO_SPMW_D_H
35 #define AIMS_IO_SPMW_D_H
36 
37 #include <cstdlib>
38 #include <aims/io/spmW.h>
39 #include <aims/io/datatypecode.h>
40 #include <aims/def/settings.h>
41 #include <aims/io/spmheader.h>
42 #include <aims/io/defaultItemW.h>
43 #include <aims/io/scaledcoding.h>
44 #include <aims/resampling/motion.h>
46 #include <stdio.h>
47 
48 namespace aims
49 {
50 
51  template <class T>
52  void SpmWriter<T>::write( const AimsData<T>& thing )
53  {
55 
56  SpmHeader hdr( thing.dimX(), thing.dimY(), thing.dimZ(), thing.dimT(),
57  thing.sizeX(), thing.sizeY(), thing.sizeZ(),
58  thing.sizeT(), _name );
59 
60  // Flip if necessary
61  const PythonHeader
62  *ph = dynamic_cast<const PythonHeader *>( thing.header() );
63  if( ph )
64  {
65  hdr.copy( *ph );
66  std::vector<int> dims(4);
67  dims[0] = thing.dimX();
68  dims[1] = thing.dimY();
69  dims[2] = thing.dimZ();
70  dims[3] = thing.dimT();
71  hdr.setProperty( "volume_dimension", dims );
72  std::vector<float> vs(4);
73  vs[0] = thing.sizeX();
74  vs[1] = thing.sizeY();
75  vs[2] = thing.sizeZ();
76  vs[3] = thing.sizeT();
77  hdr.setProperty( "voxel_size", vs );
78  // remove specific attribute that may be incompatible
79  if( hdr.hasProperty( "SPM_data_type" ) )
80  hdr.removeProperty( "SPM_data_type" );
81  if( hdr.hasProperty( "disk_data_type" ) )
82  hdr.removeProperty( "disk_data_type" );
83  if( hdr.hasProperty( "possible_data_types" ) )
84  hdr.removeProperty( "possible_data_types" );
85  }
86 
87  std::string code = carto::DataTypeCode<T>().dataType();
88  hdr.setupWriteAttributes( code );
89 
90  int y, z, x, dx, f;
91  T vmin = thing( 0 ), vmax = thing( 0 );
92  bool shen = false;
93  float offset = 0, scale = 1;
94 
95  if( code == "FLOAT" || code == "DOUBLE" )
96  {
97  double maxm = 0;
98  shen = canEncodeAsScaledS16( *thing.volume(), scale, offset, false,
99  &maxm );
100  if( shen )
101  {
102  std::cout << "16 bit coding possible: scale: " << scale
103  << ", max error: " << maxm << std::endl;
104  hdr.setProperty( "disk_data_type",
105  carto::DataTypeCode<int16_t>().dataType()
106  );
107  hdr.setProperty( "scale_factor_applied", int(0) );
108  hdr.setProperty( "scale_factor", scale );
109  }
110  else
111  std::cout << "matching interval not found. Not 16 bit codable\n";
112  }
113  ForEach4d( thing, x, y, z, f )
114  {
115  if( thing( x, y, z, f ) < vmin )
116  vmin = thing( x, y, z, f );
117  if( thing( x, y, z, f ) > vmax )
118  vmax = thing( x, y, z, f );
119  }
120 
121  if( shen )
122  hdr.setProperty( "bits_allocated", (int) 16 );
123  else
124  hdr.setProperty( "bits_allocated", (int) ( sizeof( T ) * 8 ) );
125  hdr.setProperty( "minimum", (int) vmin ) ;
126  hdr.setProperty( "maximum", (int) vmax ) ;
127 
128  const Settings sett = Settings::settings();
129  bool write4d = true;
130  try
131  {
132  write4d = (bool)
133  sett.getProperty( "spm_output_4d_volumes" )->getScalar();
134  }
135  catch( std::exception & )
136  {
137  }
138 
139  // prepare header without writing it
140  hdr.write( false, write4d );
141 
142  std::vector<float> vstom;
143  Motion stom;
144  if( hdr.getProperty( "storage_to_memory", vstom ) )
145  stom = Motion( vstom );
146  else
147  {
148  std::cout << "no storage_to_memory in SPM header\n";
149  // no storage orientation: take a default one
150  stom.rotation()(1,1) = -1;
151  stom.rotation()(2,2) = -1;
152  stom.translation()[1] = thing.dimY() - 1;
153  stom.translation()[2] = thing.dimZ() - 1;
154  bool radio = true;
155  try
156  {
157  radio = (bool) hdr.getProperty( "spm_radio_convention" )->getScalar();
158  }
159  catch( std::exception & )
160  {
161  }
162  if( !radio )
163  {
164  stom.rotation()(0,0) = -1;
165  stom.translation()[0] = thing.dimX() - 1;
166  }
167  }
168 
169  Motion m2s = stom.inverse();
170  Point3df df = m2s.transform( Point3df( thing.dimX(), thing.dimY(),
171  thing.dimZ() ) )
172  - m2s.transform( Point3df( 0, 0, 0 ) );
173  Point4d tdims = Point4d( short( rint( fabs( df[0] ) ) ),
174  short( rint( fabs( df[1] ) ) ),
175  short( rint( fabs( df[2] ) ) ), thing.dimT() );
176  Point3df incf = stom.transform( Point3df( 1, 0, 0 ) )
177  - stom.transform( Point3df( 0, 0, 0 ) );
178  Point4d inc = Point4d( int16_t( rint( incf[0] ) ),
179  int16_t( rint( incf[1] ) ),
180  int16_t( rint( incf[2] ) ), 0 );
181 
182  dx = tdims[0];
183  std::vector<int16_t> datashort;
185  if( shen )
186  datashort.insert( datashort.end(), dx, 0 );
187  std::vector<T> data( dx );
188 
189  Point3df d0f;
190  Point4d d0;
191 
192  if( write4d || thing.dimT() == 1 )
193  {
194  std::string fname = removeExtension( _name );
195  hdr.setName( fname );
196  hdr.write();
197 
198  std::string name = removeExtension( fname ) + ".img";
199  _os.open( name.c_str(), std::ios::out | std::ios::binary );
200  _os.unsetf( std::ios::skipws );
201 
202  for ( f=0; f < tdims[3]; f++)
203  for ( z=0;z<tdims[2];z++)
204  {
205  for ( y=0;y<tdims[1];y++)
206  {
207  d0f = stom.transform( Point3df( 0, y, z ) );
208  d0 = Point4d( int16_t( rint( d0f[0] ) ), int16_t( rint( d0f[1] ) ),
209  int16_t( rint( d0f[2] ) ), f );
210  if( shen )
211  {
212  for( x=0; x < dx; ++x, d0+=inc )
213  datashort[ x ] = (int16_t) rint( thing( d0 ) / scale );
214  itemWs.write( _os, &datashort[0] , dx );
215  }
216  else if( inc == Point4d( 1, 0, 0, 0 ) )
217  itemW.write( _os, &thing( d0 ) , dx );
218  else
219  {
220  for( x=0; x < dx; ++x, d0+=inc )
221  data[ x ] = thing( d0 );
222  itemW.write( _os, &data[0] , dx );
223  }
224  }
225  }
226 
227  _os.close();
228  }
229  else
230  {
231  char sequence[16];
232  int f, nt = thing.dimT();
233  std::vector<std::string> fnames;
234  std::string dname, bname;
235  dname = carto::FileUtil::dirname( _name )
237  bname = carto::FileUtil::basename( removeExtension( _name ) );
238 
239  fnames.reserve( nt );
240  for( f=0; f<nt; ++f )
241  {
242  sprintf( sequence, "%04d", f );
243  fnames.push_back( bname + std::string( sequence ) + ".img" );
244  }
245  hdr.setProperty( "series_filenames", fnames );
246 
247  // Do the job for each frame
248  for( f=0; f < nt; ++f )
249  {
250  std::string name = dname + fnames[f];
251  hdr.setName( name );
252  hdr.write( f == 0, false );
253 
254  _os.open( name.c_str(), std::ios::out | std::ios::binary );
255  _os.unsetf( std::ios::skipws );
256 
257  for (int z=0;z<tdims[2];z++)
258  for (int y=0;y<tdims[1];y++)
259  {
260  d0f = stom.transform( Point3df( 0, y, z ) );
261  d0 = Point4d( int16_t( rint( d0f[0] ) ), int16_t( rint( d0f[1] ) ),
262  int16_t( rint( d0f[2] ) ), f );
263  if( shen )
264  {
265  for( x=0; x < dx; ++x, d0+=inc )
266  datashort[ x ] = (int16_t) rint( thing( d0 ) / scale );
267  itemWs.write( _os, &datashort[0] , dx );
268  }
269  else if( inc == Point4d( 1, 0, 0, 0 ) )
270  itemW.write( _os, &thing( d0 ) , dx );
271  else
272  {
273  for( x=0; x < dx; ++x, d0+=inc )
274  data[ x ] = thing( d0 );
275  itemW.write( _os, &data[0] , dx );
276  }
277  }
278 
279  _os.close();
280  }
281  }
282  }
283 
284 }
285 
286 #endif
SPM Header class.
Definition: spmheader.h:53
Default low-levels writers.
Definition: defaultItemW.h:58
Attributed python-like header, stores all needed information about an object, currently used for volu...
Definition: pheader.h:51
int nt
int dimX() const
static std::string basename(const std::string &)
float sizeY() const
int dimT() const
virtual void copy(const PythonHeader &, bool keepUuid=false)
virtual void write(std::ostream &os, const T &item) const
Definition: defaultItemW.h:64
#define ForEach4d(thing, x, y, z, t)
static char separator()
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 &)
float sizeT() const
bool canEncodeAsScaledS16(const carto::Volume< T > &vol, float &slope, float &offset, bool enableoffset=true, double *maxerr=0)
Checks if a volume can be encoded as 16 bit signed ints with a scale factor and optionally an offset...
Definition: scaledcoding.h:69
float sizeZ() const
carto::rc_ptr< carto::Volume< T > > volume()
static Settings & settings()
float sizeX() const
void write(const AimsData< T > &thing)
Write the data with "name" file name to disk.
Definition: spmW_d.h:52
int dimZ() const
int dimY() const
const aims::Header * header() const
AffineTransformation3d inverse() const
Affine 3D transformation.