53 #include <QDataStream>
54 #include <QTextStream>
64 using namespace FSLIB;
77 , m_vecOffset(Vector3f::Zero(3))
89 , m_vecOffset(Vector3f::Zero(3))
97 Surface::Surface(
const QString &subject_id, qint32 hemi,
const QString &surf,
const QString &subjects_dir)
102 , m_vecOffset(Vector3f::Zero(3))
115 , m_vecOffset(Vector3f::Zero(3))
137 m_matTris.resize(0,3);
147 printf(
"\tcomputing normals\n");
149 MatrixX3f r1(tris.rows(),3); MatrixX3f r2(tris.rows(),3); MatrixX3f r3(tris.rows(),3);
151 for(qint32 i = 0; i < tris.rows(); ++i)
153 r1.row(i) = rr.row(
tris(i, 0));
154 r2.row(i) = rr.row(
tris(i, 1));
155 r3.row(i) = rr.row(
tris(i, 2));
158 MatrixX3f x = r2 - r1;
159 MatrixX3f y = r3 - r1;
160 MatrixX3f tri_nn(x.rows(),y.cols());
161 tri_nn.col(0) = x.col(1).cwiseProduct(y.col(2)) - x.col(2).cwiseProduct(y.col(1));
162 tri_nn.col(1) = x.col(2).cwiseProduct(y.col(0)) - x.col(0).cwiseProduct(y.col(2));
163 tri_nn.col(2) = x.col(0).cwiseProduct(y.col(1)) - x.col(1).cwiseProduct(y.col(0));
166 MatrixX3f tmp = tri_nn.cwiseProduct(tri_nn);
167 VectorXf normSize = tmp.rowwise().sum();
168 normSize = normSize.cwiseSqrt();
170 for(qint32 i = 0; i < normSize.size(); ++i)
172 tri_nn.row(i) /= normSize(i);
174 MatrixX3f
nn = MatrixX3f::Zero(rr.rows(), 3);
176 for(qint32 p = 0; p < tris.rows(); ++p)
178 Vector3i verts = tris.row(p);
179 for(qint32 j = 0; j < verts.size(); ++j)
180 nn.row(verts(j)) = tri_nn.row(p);
183 tmp = nn.cwiseProduct(nn);
184 normSize = tmp.rowwise().sum();
185 normSize = normSize.cwiseSqrt();
187 for(qint32 i = 0; i < normSize.size(); ++i)
189 nn.row(i) /= normSize(i);
197 bool Surface::read(
const QString &subject_id, qint32 hemi,
const QString &surf,
const QString &subjects_dir,
Surface &p_Surface,
bool p_bLoadCurvature)
199 if(hemi != 0 && hemi != 1)
202 QString p_sFile = QString(
"%1/%2/surf/%3.%4").arg(subjects_dir).arg(subject_id).arg(hemi == 0 ?
"lh" :
"rh").arg(surf);
204 return read(p_sFile, p_Surface, p_bLoadCurvature);
210 bool Surface::read(
const QString &path, qint32 hemi,
const QString &surf,
Surface &p_Surface,
bool p_bLoadCurvature)
212 if(hemi != 0 && hemi != 1)
215 QString p_sFile = QString(
"%1/%2.%3").arg(path).arg(hemi == 0 ?
"lh" :
"rh").arg(surf);
217 return read(p_sFile, p_Surface, p_bLoadCurvature);
227 QFile t_File(p_sFile);
229 if (!t_File.open(QIODevice::ReadOnly))
231 printf(
"\tError: Couldn't open the surface file\n");
235 printf(
"Reading surface...\n");
238 qint32 t_NameIdx = 0;
239 if(p_sFile.contains(
"lh."))
240 t_NameIdx = p_sFile.indexOf(
"lh.");
241 else if(p_sFile.contains(
"rh."))
242 t_NameIdx = p_sFile.indexOf(
"rh.");
246 p_Surface.m_sFilePath = p_sFile.mid(0,t_NameIdx);
247 p_Surface.m_sFileName = p_sFile.mid(t_NameIdx,p_sFile.size()-t_NameIdx);
249 QDataStream t_DataStream(&t_File);
250 t_DataStream.setByteOrder(QDataStream::BigEndian);
258 qint32 NEW_QUAD_FILE_MAGIC_NUMBER = 16777213;
259 qint32 TRIANGLE_FILE_MAGIC_NUMBER = 16777214;
260 qint32 QUAD_FILE_MAGIC_NUMBER = 16777215;
262 qint32 magic = IOUtils::fread3(t_DataStream);
268 if(magic == QUAD_FILE_MAGIC_NUMBER || magic == NEW_QUAD_FILE_MAGIC_NUMBER)
270 qint32 nvert = IOUtils::fread3(t_DataStream);
271 qint32 nquad = IOUtils::fread3(t_DataStream);
272 if(magic == QUAD_FILE_MAGIC_NUMBER)
273 printf(
"\t%s is a quad file (nvert = %d nquad = %d)\n", p_sFile.toLatin1().constData(),nvert,nquad);
275 printf(
"\t%s is a new quad file (nvert = %d nquad = %d)\n", p_sFile.toLatin1().constData(),nvert,nquad);
278 verts.resize(nvert, 3);
279 if(magic == QUAD_FILE_MAGIC_NUMBER)
282 for(qint32 i = 0; i < nvert; ++i)
284 for(qint32 j = 0; j < 3; ++j)
286 t_DataStream >> iVal;
287 IOUtils::swap_short(iVal);
288 verts(i,j) = ((float)iVal) / 100;
294 t_DataStream.readRawData((
char *)verts.data(), nvert*3*
sizeof(float));
295 for(qint32 i = 0; i < nvert; ++i)
296 for(qint32 j = 0; j < 3; ++j)
297 IOUtils::swap_floatp(&verts(i,j));
300 MatrixXi quads = IOUtils::fread3_many(t_DataStream, nquad*4);
301 MatrixXi quads_new(4, nquad);
303 for(qint32 j = 0; j < quads.cols(); ++j)
305 for(qint32 i = 0; i < quads.rows(); ++i)
307 quads_new(i,j) = quads(count, 0);
311 quads = quads_new.transpose();
315 faces = MatrixXi::Zero(2*nquad,3);
317 for(qint32 k = 0; k < nquad; ++k)
319 RowVectorXi quad = quads.row(k);
320 if ((quad[0] % 2) == 0)
322 faces(nface,0) = quad[0];
323 faces(nface,1) = quad[1];
324 faces(nface,2) = quad[3];
327 faces(nface,0) = quad[2];
328 faces(nface,1) = quad[3];
329 faces(nface,2) = quad[1];
334 faces(nface,0) = quad(0);
335 faces(nface,1) = quad(1);
336 faces(nface,2) = quad(2);
339 faces(nface,0) = quad(0);
340 faces(nface,1) = quad(2);
341 faces(nface,2) = quad(3);
346 else if(magic == TRIANGLE_FILE_MAGIC_NUMBER)
348 QString s = t_File.readLine();
351 t_DataStream >> nvert;
352 t_DataStream >> nface;
353 IOUtils::swap_int(nvert);
354 IOUtils::swap_int(nface);
356 printf(
"\t%s is a triangle file (nvert = %d ntri = %d)\n", p_sFile.toLatin1().constData(), nvert, nface);
357 printf(
"\t%s", s.toLatin1().constData());
360 verts.resize(3, nvert);
361 t_DataStream.readRawData((
char *)verts.data(), nvert*3*
sizeof(float));
362 for(qint32 i = 0; i < 3; ++i)
363 for(qint32 j = 0; j < nvert; ++j)
364 IOUtils::swap_floatp(&verts(i,j));
367 faces.resize(nface, 3);
369 for(qint32 i = 0; i < nface; ++i)
371 for(qint32 j = 0; j < 3; ++j)
373 t_DataStream >> iVal;
374 IOUtils::swap_int(iVal);
381 qWarning(
"Bad magic number (%d) in surface file %s",magic,p_sFile.toLatin1().constData());
385 verts.transposeInPlace();
386 verts.array() *= 0.001f;
388 p_Surface.m_matRR = verts.block(0,0,verts.rows(),3);
389 p_Surface.m_matTris = faces.block(0,0,faces.rows(),3);
392 p_Surface.m_matNN =
compute_normals(p_Surface.m_matRR, p_Surface.m_matTris);
395 if(t_File.fileName().contains(
"lh."))
396 p_Surface.m_iHemi = 0;
397 else if(t_File.fileName().contains(
"rh."))
398 p_Surface.m_iHemi = 1;
401 p_Surface.m_iHemi = -1;
406 p_Surface.m_sSurf = t_File.fileName().mid((t_NameIdx+3),t_File.fileName().size() - (t_NameIdx+3));
411 QString t_sCurvatureFile = QString(
"%1%2.curv").arg(p_Surface.m_sFilePath).arg(p_Surface.m_iHemi == 0 ?
"lh" :
"rh");
417 printf(
"\tRead a surface with %d vertices from %s\n[done]\n",nvert,p_sFile.toLatin1().constData());
425 VectorXf Surface::read_curv(
const QString &p_sFileName)
429 printf(
"Reading curvature...");
430 QFile t_File(p_sFileName);
432 if (!t_File.open(QIODevice::ReadOnly))
434 printf(
"\tError: Couldn't open the curvature file\n");
438 QDataStream t_DataStream(&t_File);
439 t_DataStream.setByteOrder(QDataStream::BigEndian);
441 qint32 vnum = IOUtils::fread3(t_DataStream);
442 qint32 NEW_VERSION_MAGIC_NUMBER = 16777215;
444 if(vnum == NEW_VERSION_MAGIC_NUMBER)
446 qint32 fnum, vals_per_vertex;
447 t_DataStream >> vnum;
449 t_DataStream >> fnum;
450 t_DataStream >> vals_per_vertex;
452 curv.resize(vnum, 1);
453 t_DataStream.readRawData((
char *)curv.data(), vnum*
sizeof(float));
454 for(qint32 i = 0; i < vnum; ++i)
455 IOUtils::swap_floatp(&
curv(i));
459 qint32 fnum = IOUtils::fread3(t_DataStream);
462 curv.resize(vnum, 1);
463 for(qint32 i = 0; i < vnum; ++i)
465 t_DataStream >> iVal;
466 IOUtils::swap_short(iVal);
467 curv(i) = ((float)iVal) / 100;
const VectorXf & curv() const
const MatrixX3f & nn() const
static MatrixX3f compute_normals(const MatrixX3f &rr, const MatrixX3i &tris)
static VectorXf read_curv(const QString &p_sFileName)
const MatrixX3i & tris() const
Surface class declaration.
static bool read(const QString &subject_id, qint32 hemi, const QString &surf, const QString &subjects_dir, Surface &p_Surface, bool p_bLoadCurvature=true)
IOUtils class declaration.