1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2008-2015 Imperial College London
5 * Copyright 2008-2015 Daniel Rueckert, Julia Schnabel
6 *
7 * Licensed under the Apache License, Version 2.0 (the "License");
8 * you may not use this file except in compliance with the License.
9 * You may obtain a copy of the License at
10 *
11 * http://www.apache.org/licenses/LICENSE-2.0
12 *
13 * Unless required by applicable law or agreed to in writing, software
14 * distributed under the License is distributed on an "AS IS" BASIS,
15 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 * See the License for the specific language governing permissions and
17 * limitations under the License.
18 */
19
20 #include "mirtk/Vector.h"
21
22 #include "mirtk/Assert.h"
23 #include "mirtk/Memory.h"
24 #include "mirtk/Array.h"
25 #include "mirtk/Indent.h"
26 #include "mirtk/Cfstream.h"
27
28 #include "mirtk/NumericsConfig.h"
29 #if MIRTK_Numerics_WITH_MATLAB
30 # include "mirtk/Matlab.h"
31 #endif
32
33
34 namespace mirtk {
35
36
37 // -----------------------------------------------------------------------------
PermuteRows(Array<int> idx)38 void Vector::PermuteRows(Array<int> idx)
39 {
40 mirtkAssert(idx.size() <= static_cast<size_t>(_rows), "valid permutation");
41 for (int r1 = 0; r1 < static_cast<int>(idx.size()); ++r1) {
42 int r2 = idx[r1];
43 if (r2 == r1) continue;
44 swap(_vector[r1], _vector[r2]);
45 for (int r = r1 + 1; r < _rows; ++r) {
46 if (idx[r] == r1) {
47 swap(idx[r], idx[r1]);
48 break;
49 }
50 }
51 }
52 }
53
54 // -----------------------------------------------------------------------------
operator <<(ostream & os,const Vector & v)55 ostream &operator <<(ostream &os, const Vector &v)
56 {
57 os << "irtkVector " << max(0, v._rows) << endl;
58 if (v._rows > 0) {
59 const bool swapped = (GetByteOrder() == LittleEndian);
60 if (swapped) {
61 swap64((char *)v._vector, (char *)v._vector, v._rows);
62 }
63 os.write((char *)v._vector, v._rows * sizeof(double));
64 if (swapped) {
65 swap64((char *)v._vector, (char *)v._vector, v._rows);
66 }
67 }
68 return os;
69 }
70
71 // -----------------------------------------------------------------------------
operator >>(istream & is,Vector & v)72 istream &operator >> (istream &is, Vector &v)
73 {
74 // Read file type identifier
75 char keyword[11];
76 is.read(keyword, 11);
77 if (strncmp(keyword, "irtkVector ", 11) != 0) {
78 cerr << "Vector: File does not appear to be a binary (M)IRTK Vector file" << endl;
79 exit(1);
80 }
81
82 // Read vector size
83 int rows = 0;
84 is >> rows;
85 if (!is || rows < 0) {
86 cerr << "Vector: Failed to read vector size from binary file" << endl;
87 exit(1);
88 }
89
90 // Skip remaining characters (e.g., comment) on first line
91 char buffer[256];
92 while (is.get() != '\n') {
93 is.get(buffer, 256);
94 }
95
96 // Read vector data
97 v.Initialize(rows);
98 if (rows > 0) {
99 is.read((char *)v._vector, rows * sizeof(double));
100 if (is.fail()) {
101 cerr << "Vector: File contains fewer values than expected: #rows = " << rows << endl;
102 exit(1);
103 }
104 if (GetByteOrder() == LittleEndian) {
105 swap64((char *)v._vector, (char *)v._vector, v._rows);
106 }
107 }
108
109 return is;
110 }
111
112 // -----------------------------------------------------------------------------
operator <<(Cofstream & to,const Vector & v)113 Cofstream &operator <<(Cofstream& to, const Vector &v)
114 {
115 to.WriteAsChar("irtkVector", 11);
116 const int rows = max(0, v._rows);
117 to.WriteAsInt(&rows, 1);
118 if (v._rows > 0) to.WriteAsDouble(v._vector, v._rows);
119 return to;
120 }
121
122 // -----------------------------------------------------------------------------
operator >>(Cifstream & from,Vector & v)123 Cifstream &operator >>(Cifstream &from, Vector &v)
124 {
125 char keyword[11];
126 from.ReadAsChar(keyword, 11);
127 keyword[10] = '\0';
128 if (strcmp(keyword, "irtkVector") != 0) {
129 cerr << "Vector: File does not appear to be a binary (M)IRTK Vector file" << endl;
130 exit(1);
131 }
132
133 int rows = 0;
134 if (!from.ReadAsInt(&rows, 1) || rows < 0) {
135 cerr << "Vector: Failed to read vector size from binary file" << endl;
136 exit(1);
137 }
138
139 v.Initialize(rows);
140 if (rows > 0 && !from.ReadAsDouble(v._vector, rows)) {
141 cerr << "Vector: File contains fewer values than expected: #rows = " << rows << endl;
142 exit(1);
143 }
144
145 return from;
146 }
147
148 // -----------------------------------------------------------------------------
Print(Indent indent) const149 void Vector::Print(Indent indent) const
150 {
151 cout << indent << "Vector " << _rows << endl;
152 ++indent;
153 cout.setf(ios::right);
154 cout.setf(ios::fixed);
155 cout.precision(4);
156 for (int i = 0; i < _rows; i++) {
157 cout << indent << setw(15) << _vector[i] << endl;
158 }
159 cout.precision(6);
160 cout.unsetf(ios::right);
161 cout.unsetf(ios::fixed);
162 }
163
164 // -----------------------------------------------------------------------------
Read(const char * filename)165 void Vector::Read(const char *filename)
166 {
167 // Open file stream
168 ifstream from(filename, ios::in | ios::binary);
169
170 // Check whether file opened ok
171 if (!from) {
172 cerr << "Vector::Read: Can't open file " << filename << endl;
173 exit(1);
174 }
175
176 // Read vector
177 from >> *this;
178 }
179
180 // -----------------------------------------------------------------------------
Write(const char * filename) const181 void Vector::Write(const char *filename) const
182 {
183 // Open file stream
184 ofstream to(filename, ios::out | ios::binary);
185
186 // Check whether file opened ok
187 if (!to) {
188 cerr << "Vector::Write: Can't open file " << filename << endl;
189 exit(1);
190 }
191
192 // Write vector
193 to << *this;
194 }
195
196 #if MIRTK_Numerics_WITH_MATLAB
197
198 // -----------------------------------------------------------------------------
VectorToMxArray(const Vector & v)199 inline mxArray *VectorToMxArray(const Vector &v)
200 {
201 Matlab::Initialize();
202 mxArray *m = mxCreateDoubleMatrix(v.Rows(), 1, mxREAL);
203 memcpy(mxGetPr(m), v.RawPointer(), v.Rows() * sizeof(double));
204 return m;
205 }
206
207 // -----------------------------------------------------------------------------
WriteMAT(const char * fname,const char * varname) const208 bool Vector::WriteMAT(const char *fname, const char *varname) const
209 {
210 Matlab::Initialize();
211 MATFile *fp = matOpen(fname, "w");
212 if (fp == NULL) return false;
213 mxArray *m = VectorToMxArray(*this);
214 if (matPutVariable(fp, varname, m) != 0) {
215 mxDestroyArray(m);
216 matClose(fp);
217 return false;
218 }
219 mxDestroyArray(m);
220 return (matClose(fp) == 0);
221 }
222
223 #endif // MIRTK_Numerics_WITH_MATLAB
224
225
226 } // namespace mirtk
227