8

I want to read the not reconstructed data from the Stanford Bunny. The point data is stored as several range images, which have to be transformed to be combined to one big point cloud, like written in the README:

These data files were obtained with a Cyberware 3030MS optical
triangulation scanner.  They are stored as range images in the "ply"
format.  The ".conf" file contains the transformations required to
bring each range image into a single coordinate system.

This is the .conf-file:

camera -0.0172 -0.0936 -0.734  -0.0461723 0.970603 -0.235889 0.0124573
bmesh bun000.ply 0 0 0 0 0 0 1
bmesh bun045.ply -0.0520211 -0.000383981 -0.0109223 0.00548449 -0.294635 -0.0038555 0.955586
bmesh bun090.ply 2.20761e-05 -3.34606e-05 -7.20881e-05 0.000335889 -0.708202 0.000602459 0.706009
bmesh bun180.ply 0.000116991 2.47732e-05 -4.6283e-05 -0.00215148 0.999996 -0.0015001 0.000892527
bmesh bun270.ply 0.000130273 1.58623e-05 0.000406764 0.000462632 0.707006 -0.00333301 0.7072
bmesh top2.ply -0.0530127 0.138516 0.0990356 0.908911 -0.0569874 0.154429 0.383126
bmesh top3.ply -0.0277373 0.0583887 -0.0796939 0.0598923 0.670467 0.68082 -0.28874
bmesh bun315.ply -0.00646017 -1.36122e-05 -0.0129064 0.00449209 0.38422 -0.00976512 0.923179
bmesh chin.ply 0.00435102 0.0882863 -0.108853 -0.441019 0.213083 0.00705734 0.871807
bmesh ear_back.ply -0.0829384 0.0353082 0.0711536 0.111743 0.925689          -0.215443 -0.290169

For each range image seven values are stored. But I do not know, what information can be obtained from these values. I guess that three of them will contain some information about the translation and maybe three contain information about the rotation. But I didn't find something about the order of these values and how to transform the values to get one point cloud.

The wiki page doesn't handle with range images and I found nothing more at the Stanford pages. They just talk about, that the method of Turk94 is used to scan this data set, but the method has no information about the transformations needed. (Or I was not able to get the information out of this paper.)

Does anybody know how to read these values correctly? Why is there a transformation for the camera position? Is this just a good initial value to view the whole point cloud?

Thanks for your help.

EDIT:

Ok. At this point, I already tried to read the data and to correctly transform them, but everything did not work. I use the boost library to handle with the quaternions

Here is my code for it:

boost::math::quaternion<double> translation, quaternionRotation;
//Get Transformation
translation = boost::math::quaternion<double>(0.0, lineData[2].toDouble(), lineData[3].toDouble(), lineData[4].toDouble());
quaternionRotation = boost::math::quaternion<double>(lineData[5].toDouble(),lineData[6].toDouble(),lineData[7].toDouble(),lineData[8].toDouble());

//do some file related stuff 
//...

//for each line: read the point data and transform it and store the point in a data array
pointData[j].x = stringPointData[0].toDouble();
pointData[j].y = stringPointData[1].toDouble();
pointData[j].z = stringPointData[2].toDouble();
tmpQuat = boost::math::quaternion<double> (0.0,pointData[j].x,pointData[j].y,pointData[j].z);
//first translation
tmpQuat += translation;
//then quaternion rotation
tmpQuat = (quaternionRotation * (tmpQuat) * boost::math::conj(quaternionRotation));
//read the data from quaternion to a usual type
pointData[j].x = tmpQuat.R_component_2();
pointData[j].y = tmpQuat.R_component_3();
pointData[j].z = tmpQuat.R_component_4();

I assume that the first component of the quaternion is the w component and the others refers to x, y andz like in equation 2 from here. If necessary I can provide the screenshots of the false transformations.

EDIT: It is written in the source code of zipper in the file zipper.c, that the 7 values are saved as followed:

transX transY transZ quatX quatY quatZ quatW

The quaternion is then transformed into a rotation matrix and then the rotation is performed with this new matrix. But even with this information, I am not able to transform it correctly. To test it, I implemented the function quat_to_mat() from zipper in my project:

glm::dmat4 cPlyObjectLoader::quat_to_mat(boost::math::quaternion<double> quat) const
{
 float s;
 float xs,ys,zs;
 float wx,wy,wz;
 float xx,xy,xz;
 float yy,yz,zz;
 glm::dmat4 mat(1.0);

 s = 2 / (quat.R_component_2()*quat.R_component_2() +
          quat.R_component_3()*quat.R_component_3() + 
          quat.R_component_4()*quat.R_component_4() + 
          quat.R_component_1()*quat.R_component_1());

 xs = quat.R_component_2() * s;
 ys = quat.R_component_3() * s;
 zs = quat.R_component_4() * s;

 wx = quat.R_component_1() * xs;
 wy = quat.R_component_1() * ys;
 wz = quat.R_component_1() * zs;

 xx = quat.R_component_2() * xs;
 xy = quat.R_component_2() * ys;
 xz = quat.R_component_2() * zs; 

 yy = quat.R_component_3() * ys;
 yz = quat.R_component_3() * zs;
 zz = quat.R_component_4() * zs;

 mat[0][0] = 1 - (yy + zz);
 mat[0][1] = xy - wz;
 mat[0][2] = xz + wy;
 mat[0][3] = 0;

 mat[1][0] = xy + wz;
 mat[1][1] = 1 - (xx + zz);
 mat[1][2] = yz - wx;
 mat[1][3] = 0;

 mat[2][0] = xz - wy;
 mat[2][1] = yz + wx;
 mat[2][2] = 1 - (xx + yy);
 mat[2][3] = 0;

 mat[3][0] = 0;
 mat[3][1] = 0;
 mat[3][2] = 0;
 mat[3][3] = 1;

 return mat;
}

Now I am doing the translation and rotation with a vector and this matrix:

quaternionRotation = boost::math::quaternion<double>(lineData[8].toDouble(),lineData[5].toDouble(),lineData[6].toDouble(),lineData[7].toDouble());
rotationMat = this->quat_to_mat(quaternionRotation);
translationVec = glm::dvec4(lineData[2].toDouble(), lineData[3].toDouble(), lineData[4].toDouble(),0.0);

//same stuff as above
//...

glm::dvec4 curPoint =   glm::dvec4(pointData[j].x,pointData[j].y,pointData[j].z,1.0);
curPoint += translationVec;
curPoint = rotationMat*curPoint;

The result is different to my quaternion rotation (Why? It should be the same.), but not correct.

Debug information:

  • the input of all transformations is correct
  • the input of all points is correct
1
  • Ok. It seems, that the first three values refers to a (x,y,z)-translation and the last four values refer to a quaternion rotation. If nobody else knows something about it, I will write a whole answer after some test if this is really correct.
    – DanceIgel
    Jun 11, 2015 at 15:15

6 Answers 6

7

As i read from stanford 3d scan

For all the Stanford models, alignment was done using a modified ICP algorithm, as described in this paper. These alignments are stored in ".conf" files, which list each range image in the model along with a translation and a quaternion rotation.

Here is the link to "this paper"

Edit: The two methods are called zippering and volmetric merging

6
  • This paper is exactly the one from my question (Turk94) and I am not able to get the transformation information out of this paper.
    – DanceIgel
    Jun 11, 2015 at 14:37
  • Oh my bad, though it stands in the pdf file. Did you read the volumetric merging paper as well?
    – Ello
    Jun 11, 2015 at 14:38
  • Yes, that's correct, but If we assume that the first three values refer to a translation and the four other values refer to the quaternion rotation, I am still not able to handle with quaternion rotation. But I will read about it and will write new stuff later, after some tests. Thanks for your answer!
    – DanceIgel
    Jun 11, 2015 at 14:47
  • I did not read the volumetric merging paper and this should not be necessary, cause the transformation values are given in the .conf-file. But it still does not work, even with a working quaternion rotation.
    – DanceIgel
    Jul 21, 2015 at 15:29
  • I am sorry that I can not help you any further with this. I will ask a friend of mine, maybe he has some information
    – Ello
    Jul 21, 2015 at 15:33
7

As Ello mentioned, it is written at the stanford 3D repo:

For all the Stanford models, alignment was done using a modified ICP algorithm, as described in this paper. These alignments are stored in ".conf" files, which list each range image in the model along with a translation and a quaternion rotation.

But that is not enough to understand everything of this data file.

It is correct, that the first line:

camera -0.0172 -0.0936 -0.734  -0.0461723 0.970603 -0.235889 0.0124573

stores a good initial camera position and every other line starting with bmesh refers to a .ply-file, which stores a ranged image.

The transformation values are stored as followed:

transX transY transZ quatX quatY quatZ quatW

where trans... refers to a translation value and quat... refers to a value of the quaternion. Currently, I do not know, why it doesn't work with the quaternion rotation by itself, but by transforming it into a rotation matrix with the code of zipper the transformation is correct. Be aware, that the translation is stored first, but to get a correct transformation the rotation has to be done at the beginning and the translation afterwards.

My code snippet to read the files and transform it, is the following:

boost::math::quaternion<double> translation, quaternionRotation;
//Get Transformation
translationVec = glm::dvec4(lineData[2].toDouble(), lineData[3].toDouble(), lineData[4].toDouble(),0.0);
quaternionRotation = boost::math::quaternion<double>(lineData[8].toDouble(),lineData[5].toDouble(),lineData[6].toDouble(),lineData[7].toDouble());
//calculate the unit quaternion
double magnitude = std::sqrt(
      quaternionRotation.R_component_1()*quaternionRotation.R_component_1()+
      quaternionRotation.R_component_2()*quaternionRotation.R_component_2()+
      quaternionRotation.R_component_3()*quaternionRotation.R_component_3()+
      quaternionRotation.R_component_4()*quaternionRotation.R_component_4());
quaternionRotation /= magnitude;
rotationMat = this->quat_to_mat(quaternionRotation);

//do some file related stuff 
//...

//for each line: read the point data and transform it and store the point in a data array
pointData[j].x = stringPointData[0].toDouble();
pointData[j].y = stringPointData[1].toDouble();
pointData[j].z = stringPointData[2].toDouble();
//transform the curren point
glm::dvec4 curPoint =  glm::dvec4(pointData[j].x,pointData[j].y,pointData[j].z,1.0);
//first rotation
curPoint = rotationMat*curPoint;
//then translation
curPoint += translationVec;
//store the data in a data array
pointData[j].x = curPoint.x;
pointData[j].y = curPoint.y;
pointData[j].z = curPoint.z;

I know, that it's not the best one, but it works. Feel free to optimize it by yourself.

3
  • I am glad, that you found a solution.
    – Ello
    Jul 28, 2015 at 9:03
  • I'm using Matlab (rather than PCL. It provides a function quat2angle( ... ) which converts a quaternion to a set of three rotational angles. The key here is that there are many different ways to perform the three rotations that make up the quaternion. Matlab docs list these as 'ZYX', 'ZYZ', 'ZXY', 'ZXZ', 'YXZ', 'YXY', 'YZX', 'YZY', 'XYZ', 'XYX', 'XZY', and 'XZX'. By a process of trial and error it seems that the quaternions in the conf file represent rotations sequence 'YXZ'. You would need to get the rotation sequence correct for the transformation to work. Could this be the problem ? Dec 10, 2015 at 22:24
  • +Dave Durbin No this is not the problem. It works fine with the transformation of the quaternion to a rotation matrix but it does not work with the quaternion itself.
    – DanceIgel
    Dec 11, 2015 at 13:25
0

Here is the file converter that I wrote. It will assemble all the scans into a single file, one point per line. It supports different file formats (including Stanford .conf files).

#include <string>
#include <vector>
#include <sstream>
#include <iostream>
#include <stdio.h>
#include <ctype.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#ifndef M_PI
#define M_PI 3.14159265
#endif

class LineInput {
  public:
    LineInput(const std::string& filename) {
        F_ = fopen(filename.c_str(), "r" ) ;
        ok_ = (F_ != 0) ;
    }
    ~LineInput() {
        if(F_ != 0) {
            fclose(F_); F_ = 0 ; 
        }
    }
    bool OK() const { return ok_ ; }
    bool eof() const { return feof(F_) ; }
    bool get_line() {
        line_[0] = '\0' ;
        // Skip the empty lines
        while(!isprint(line_[0])) {
            if(fgets(line_, MAX_LINE_LEN, F_) == 0) {
                return false ;
            }
        }
        // If the line ends with a backslash, append
        // the next line to the current line.
        bool check_multiline = true ;
        int total_length = MAX_LINE_LEN ;
        char* ptr = line_ ;
        while(check_multiline) {
            int L = strlen(ptr) ;
            total_length -= L ;
            ptr = ptr + L - 2;
            if(*ptr == '\\' && total_length > 0) {
                *ptr = ' ' ;
                ptr++ ;
                fgets(ptr, total_length, F_) ;
            } else {
                check_multiline = false ;
            }
        }
        if(total_length < 0) {
            std::cerr
                << "MultiLine longer than " 
                << MAX_LINE_LEN << " bytes" << std::endl ;
        }
        return true ;
    }
    int nb_fields() const { return field_.size() ;  }
    char* field(int i) { return field_[i] ;   }
    int field_as_int(int i) {
        int result ;
        ok_ = ok_ && (sscanf(field(i), "%d", &result) == 1) ;
        return result ;
    }
    double field_as_double(int i) {
        double result ;
        ok_ = ok_ && (sscanf(field(i), "%lf", &result) == 1) ;
        return result ;
    }
    bool field_matches(int i, const char* s) {
        return !strcmp(field(i), s) ;
    }
    void get_fields(const char* separators=" \t\r\n") {
        field_.resize(0) ;
        char* tok = strtok(line_,separators) ;
        while(tok != 0) { 
            field_.push_back(tok) ; 
            tok = strtok(0,separators) ;
        }
    }

  private:
    enum { MAX_LINE_LEN = 65535 } ;
    FILE* F_ ;
    char line_[MAX_LINE_LEN] ;
    std::vector<char*> field_ ;
    bool ok_ ;
} ;

std::string to_string(int x, int mindigits) {
    char buff[100] ;
    sprintf(buff, "%03d", x) ;
    return std::string(buff) ;
}

double M[4][4] ;

void transform(double* xyz) {
    double xyzw[4] ;
    for(unsigned int c=0; c<4; c++) {
        xyzw[c] = M[3][c] ;
    }
    for(unsigned int j=0; j<4; j++) {
        for(unsigned int i=0; i<3; i++) {
            xyzw[j] += M[i][j] * xyz[i] ;
        }
    }
    for(unsigned int c=0; c<3; c++) {
        xyz[c] = xyzw[c] / xyzw[3] ;
    }
}

bool read_frames_file(int no) {
    std::string filename = "scan" + to_string(no,3) + ".frames" ;
    std::cerr << "Reading frames from:" << filename << std::endl ;
    LineInput in(filename) ;
    if(!in.OK()) {
       std::cerr << " ... not found" << std::endl ;
       return false ;
    }
    while(!in.eof() && in.get_line()) {
        in.get_fields() ;
        if(in.nb_fields() == 17) {
            int f = 0 ;
            for(unsigned int i=0; i<4; i++) {
                for(unsigned int j=0; j<4; j++) {
                    M[i][j] = in.field_as_double(f) ; f++ ;
                }
            }
        } 
    }
    return true ;
}

bool read_pose_file(int no) {
    std::string filename = "scan" + to_string(no,3) + ".pose" ;
    std::cerr << "Reading pose from:" << filename << std::endl ;
    LineInput in(filename) ;
    if(!in.OK()) {
       std::cerr << " ... not found" << std::endl ;
       return false ;
    }
    double xyz[3] ;
    double euler[3] ;
    in.get_line() ;
    in.get_fields() ;
    xyz[0] = in.field_as_double(0) ;
    xyz[1] = in.field_as_double(1) ;
    xyz[2] = in.field_as_double(2) ;
    in.get_line() ;
    in.get_fields() ;
    euler[0] = in.field_as_double(0) * M_PI / 180.0 ;
    euler[1] = in.field_as_double(1) * M_PI / 180.0 ;
    euler[2] = in.field_as_double(2) * M_PI / 180.0 ;

   double sx = sin(euler[0]);
   double cx = cos(euler[0]);
   double sy = sin(euler[1]);
   double cy = cos(euler[1]);
   double sz = sin(euler[2]);
   double cz = cos(euler[2]);

   M[0][0] = cy*cz;
   M[0][1] = sx*sy*cz + cx*sz;
   M[0][2] = -cx*sy*cz + sx*sz;
   M[0][3] = 0.0;
   M[1][0] = -cy*sz;
   M[1][1] = -sx*sy*sz + cx*cz;
   M[1][2] = cx*sy*sz + sx*cz;
   M[1][3] = 0.0;
   M[2][0] = sy;
   M[2][1] = -sx*cy;
   M[2][2] = cx*cy;
   M[2][3] = 0.0;
   M[3][0] = xyz[0];
   M[3][1] = xyz[1];
   M[3][2] = xyz[2];
   M[3][3] = 1.0;
   return true ;
}

void setup_transform_from_translation_and_quaternion(
    double Tx, double Ty, double Tz,
    double Qx, double Qy, double Qz, double Qw
) {
    /* for unit q, just set s = 2 or set xs = Qx + Qx, etc. */

    double s = 2.0 / (Qx*Qx + Qy*Qy + Qz*Qz + Qw*Qw);

    double xs = Qx * s;
    double ys = Qy * s;
    double zs = Qz * s;

    double wx = Qw * xs;
    double wy = Qw * ys;
    double wz = Qw * zs;

    double xx = Qx * xs;
    double xy = Qx * ys;
    double xz = Qx * zs;

    double yy = Qy * ys;
    double yz = Qy * zs;
    double zz = Qz * zs;

    M[0][0] = 1.0 - (yy + zz);
    M[0][1] = xy - wz;
    M[0][2] = xz + wy;
    M[0][3] = 0.0;

    M[1][0] = xy + wz;
    M[1][1] = 1 - (xx + zz);
    M[1][2] = yz - wx;
    M[1][3] = 0.0;

    M[2][0] = xz - wy;
    M[2][1] = yz + wx;
    M[2][2] = 1 - (xx + yy);
    M[2][3] = 0.0;

    M[3][0] = Tx;
    M[3][1] = Ty;
    M[3][2] = Tz;
    M[3][3] = 1.0;
}

bool read_points_file(int no) {
    std::string filename = "scan" + to_string(no,3) + ".3d" ;
    std::cerr << "Reading points from:" << filename << std::endl ;
    LineInput in(filename) ;
    if(!in.OK()) {
       std::cerr << " ... not found" << std::endl ;
       return false ;
    }
    while(!in.eof() && in.get_line()) {
        in.get_fields() ;
        double xyz[3] ;
        if(in.nb_fields() >= 3) {
            for(unsigned int c=0; c<3; c++) {
                xyz[c] = in.field_as_double(c) ;
            }
            transform(xyz) ;
            printf("%f %f %f\n",xyz[0],xyz[1],xyz[2]) ;
        }
    }
    return true ;
}


/* only works for ASCII PLY files */
void read_ply_file(char* filename) {
    std::cerr << "Reading points from:" << filename << std::endl;
    LineInput in(filename) ;
    if(!in.OK()) {
        std::cerr << filename << ": could not open" << std::endl ;
        return;
    }
    bool reading_vertices = false;
    int nb_vertices = 0 ;
    int nb_read_vertices = 0 ;
    while(!in.eof() && in.get_line()) {
        in.get_fields();
        if(reading_vertices) {
            double xyz[3] ;
            for(unsigned int c=0; c<3; c++) {
                xyz[c] = in.field_as_double(c) ;
            }
            transform(xyz) ;
            printf("%f %f %f\n",xyz[0],xyz[1],xyz[2]) ;
            ++nb_read_vertices;
            if(nb_read_vertices == nb_vertices) {
                return;
            }
        } else if(
            in.field_matches(0,"element") &&
            in.field_matches(1,"vertex") 
        ) {
            nb_vertices = in.field_as_int(2);
        } else if(in.field_matches(0,"end_header")) {
            reading_vertices = true;
        }
    }    
}

/* For Stanford scanning repository */
void read_conf_file(char* filename) {
    LineInput in(filename) ;
    if(!in.OK()) {
        std::cerr << filename << ": could not open" << std::endl ;
        return;
    }
    while(!in.eof() && in.get_line()) {
        in.get_fields();
        if(in.nb_fields() == 0) { continue ; }
        if(in.field_matches(0,"bmesh")) {
            char* filename = in.field(1);
            // Translation vector
            double Tx = in.field_as_double(2);
            double Ty = in.field_as_double(3);
            double Tz = in.field_as_double(4);
            /// Quaternion
            double Qx = in.field_as_double(5);
            double Qy = in.field_as_double(6);
            double Qz = in.field_as_double(7);
            double Qw = in.field_as_double(8);
            setup_transform_from_translation_and_quaternion(Tx,Ty,Tz,Qx,Qy,Qz,Qw);
            read_ply_file(filename);
        }
    }
 }


int main(int argc, char** argv) {
    if(argc != 2) { return -1 ; }
    if(strstr(argv[1],".conf")) {
        read_conf_file(argv[1]);
    } else {
        int max_i = atoi(argv[1]) ;
        for(int i=0; i<=max_i; i++) {
            if(!read_frames_file(i)) {
                read_pose_file(i) ;
            }
            read_points_file(i) ;
        }
    }
    return 0 ;
}
10
  • Thanks for your answer! You are using the transformation from the quaternion to a rotaion matrix too. Did you tried to transform it with the orginal quaternion?
    – DanceIgel
    Jul 28, 2015 at 9:51
  • I preferred to create a matrix from the quaternion, so that my program works with different file formats (that store calibration data in matrices), but clearly it is possible to directly use the quaternion.
    – BrunoLevy
    Jul 28, 2015 at 9:53
  • Yeah, but I tried it, and it failed. Do you have any ideas, why this happens? The code should be posted above.
    – DanceIgel
    Jul 28, 2015 at 9:57
  • I'm unsure, but by quickly looking at your code, I did not see the divide by the 'w' coordinate, right after the matrix-vector product (in my code, it is in the transform() function). Tell me (I'll try digging deeper if it's not the problem).
    – BrunoLevy
    Jul 28, 2015 at 10:31
  • I have no problem with the matrix-vector rotation. The pure quaternion rotation does not work. I took the information from gamasutra.com/view/feature/131686/… to do my quaternion roation.
    – DanceIgel
    Jul 28, 2015 at 12:04
0

Okay so here is my solution since none of the above worked for me (note this is in python using blender's bpy). It seems that I need to transpose the rotation part of my 4x4 transformation matrix (note I am using a standard way to convert quaternion to rotation matrix and not the one from zipper). Also note since I am using blender when importing or using any model it only stores the models local coordinates relative to the objects world transformation so you do not have to do this point = objWorld * point, it is blender specific.

#loop
for meshName, transform in zip(plyFile, transformations):
    #Build Quaternion
    #transform structure [x, y, z, qx, qy, qz, qw]
    Rt = mathutils.Quaternion((transform[6], transform[3], transform[4], transform[5])).to_matrix().to_4x4()
    Rt.normalize()
    Rt.transpose()
    Rt[0][3] = transform[0]
    Rt[1][3] = transform[1]
    Rt[2][3] = transform[2]

    bpy.ops.object.select_all(action='DESELECT')
    #import the ply mesh into blender
    bpy.ops.import_mesh.ply(filepath=baseDir + meshName)
    #get the ply object
    obj = bpy.context.object
    #get objects world matrix
    objWorld = obj.matrix_world

    for index in range(len(obj.data.vertices)):
        #get local point
        point = mathutils.Vector([obj.data.vertices[index].co[0],obj.data.vertices[index].co[1], obj.data.vertices[index].co[2], 1.])
        #convert local point to world
        point = objWorld * point
        #apply ply transformation
        point = Rt * point
        #update the point in the mesh
        obj.data.vertices[index].co[0] = point[0]
        obj.data.vertices[index].co[1] = point[1]
        obj.data.vertices[index].co[2] = point[2]
#all vertex positions should be updated correctly
0

As mentioned in other answers, the Stanford 3D repository gives some info about the data organization in the '.conf' files but, the transformation for the bunny model were not working properly when using the quaternion data provided.

I was also stuck in this registration problem for the bunny model, and based on my tests I have some extra considerations to add up. When applying the transformation - rotations to be more specific - I have realized that quaternion values were not rotating the cloud in the correct direction but, when using the corresponding Euler notation, by changing the sign of one specific axis of rotation, I got the correct registration. So, back to the quaternion notation used in the '.conf' file, after some tests I have noticed that by changing the sign of the 'w' component in the quaternion, in each 'bmesh' row, but the first (bun000.ply), the rotation by quaternion can be used.

Furthermore, for some reason, when registering the dragon (dragon_stand and dragon_side) and armadillo (armadillo_stand) stanford point clouds, in order to get the correct result I had to use a different sequence for reading the quaternion data in the ‘.conf’ file. It seems to be stored as: tx ty tz qw qx qy qz where 't' refers to a translation value and 'q' refers to a quaternion value. Just to be clear, I have just tested these three models, therefore, I don’t know what is the default pattern for the quaternion values. Besides, for these last two point cloud models, I did not need to change the '.conf' file.

I hope this could be useful for someone else trying to do the same

0

Just in case someone is looking for a full python implementation on the basis of what @DanceIgel found out, here is some code in python 3.9.1, also generating a figure in mathplotlib:

# Python 3.9.1

import numpy as np
import sys
import math
import glob

import matplotlib  
matplotlib.use('TkAgg')   
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import open3d as o3d

def get_pointcloud_files(path):
    files = list()
    for f in glob.glob(path + '/*.ply'):
        files.append(f)
    return files

def get_pointcloud_from_file(path, filename):
    cloud = o3d.io.read_point_cloud(path + '/' + filename)
    return cloud

def get_transformations_from_file(path, filename):
    with open(path + '/' + filename) as f:
        lines = (line for line in f)
        source = np.loadtxt(lines, delimiter=' ', skiprows=1, dtype='str')
        source = np.delete(source, 0, 1)  #remove camera
        
        filenames = source[:,0]
        source = source[filenames.argsort()]
        filenames = np.sort(filenames)
        
        translations = list()
        for row in source[:,1:4]:
            translations.append(np.reshape(row, [3,1]).astype(np.float32))
        
        quaternions = list()
        for row in source[:,4:]:
            quaternions.append(np.reshape(row, [4,1]).astype(np.float32))    
        
    return filenames, translations, quaternions

def quaternion_rotation_matrix(Q):

    # Extract the values from Q
    q0 = Q[3]
    q1 = Q[0]
    q2 = Q[1]
    q3 = Q[2]
    
    # calculate unit quarternion
    magnitude = math.sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3)
    
    q0 = q0 / magnitude
    q1 = q1 / magnitude
    q2 = q2 / magnitude
    q3 = q3 / magnitude
     
    # First row of the rotation matrix
    r00 = 2 * (q0 * q0 + q1 * q1) - 1
    r01 = 2 * (q1 * q2 - q0 * q3)
    r02 = 2 * (q1 * q3 + q0 * q2)
     
    # Second row of the rotation matrix
    r10 = 2 * (q1 * q2 + q0 * q3)
    r11 = 2 * (q0 * q0 + q2 * q2) - 1
    r12 = 2 * (q2 * q3 - q0 * q1)
     
    # Third row of the rotation matrix
    r20 = 2 * (q1 * q3 - q0 * q2)
    r21 = 2 * (q2 * q3 + q0 * q1)
    r22 = 2 * (q0 * q0 + q3 * q3) - 1
     
    # 3x3 rotation matrix
    rot_matrix = np.array([[r00, r01, r02],
                           [r10, r11, r12],
                           [r20, r21, r22]])
    
    rot_matrix = np.transpose(rot_matrix)
                            
    return rot_matrix

if __name__=="__main__": # $python visualization_bunny.py bunny/data

    path = sys.argv[1]
    
    # load transformations and filenames from file
    filenames, translations, quaternions = get_transformations_from_file(path, 'bun.conf')
    
    curr_transformation = np.zeros([3,4])

    clouds = list()
    for curr_filename, curr_quaternion, curr_translation in zip(filenames, quaternions, translations):  # go through input files
        curr_cloud = get_pointcloud_from_file(path, curr_filename)
        
        # convert cloud to numpy
        curr_cloud = np.asarray(curr_cloud.points)
        
        # compute rotation matrix from quaternions
        curr_rotation_matr = quaternion_rotation_matrix(curr_quaternion)
        curr_rotation_matr = np.squeeze(curr_rotation_matr)
        curr_translation = np.squeeze(curr_translation)
        
        # create transformation matrix
        curr_transformation[:,0:3] = curr_rotation_matr
        curr_transformation[:,3] = curr_translation
        
        # transform current cloud
        for i in range(curr_cloud.shape[0]):
            # apply rotation
            curr_point = np.matmul(curr_rotation_matr, np.transpose(curr_cloud[i,:]))
            # apply translation
            curr_point = curr_point + curr_translation 
            
            curr_cloud[i,0] = curr_point[0]
            curr_cloud[i,1] = curr_point[1]
            curr_cloud[i,2] = curr_point[2]
            
        # add current cloud to list of clouds
        clouds.append(curr_cloud)
    
    #plot separate point clouds in same graph
    ax = plt.axes(projection='3d')
    for cloud in clouds:
        ax.plot(cloud[:,0], cloud[:,1], cloud[:,2], 'bo', markersize=0.005)
    #ax.view_init(elev=90, azim=270)
    ax.view_init(elev=100, azim=270)
    plt.axis('off')
    plt.savefig("ZZZ_Stanford_Bunny_PointCloud.png", bbox_inches='tight')
    plt.show()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.