Lx=d: SundayResearch

The floating point reproducibility problem On a typical computer, the results of floating point operations can differ from the mathematical result. In general, floating point operations like addition may not satisfy all mathematical properties for that operation. For example, the associativity property of addition may not hold for all numerical values: $(x + y) + z \stackrel{?}{=} x + (y + z)$ A decimal representation of a number cannot always be represented by a finite binary representation in computer hardware. For example, the computed result of $0.1 + 0.01 - 0.1 - 0.01$ is not equal to $0$, but mathematically these values should be equivalent. For example, you can demonstrate this using the Python interpreter (Python 3.8.5): >>> 0.1 + 0.01 - 0.1 - 0.01 -5.204170427930421e-18 But mathematical properties do hold in other cases. For example, the result of $0.1 - 0.1 + 0.01 - 0.01$ is equal to $0$. >>> 0.1 - 0.1 + 0.01 - 0.01 0.0 Another problem results from rounding...

Get camera parameters from the OpenGL modelview matrix OpenGL modelview matrix includes model related operations, therefore, the scene was scaled, translated, or rotated, these effects are also in there. There is no way to separate these operation from the pure camera parameters from the matrix only. Here we assume the model operation is identity. If your renderer separated the model operations and camera parameters, you need to undo the model operation with your implementation dependent way. Our renderer has separation of them, so we needed to undo the effect. The elements of OpenGL modelview matrix is the following. If you want to know the details of the each element, please refer my other blog (http://shitohichiumaya.blogspot.de/2011/01/what-matrix-glulookat-generates-1.html) \begin{eqnarray*}  \left[   \begin{array}{cccc}    x_x & y_x & z_x & 0 \\    x_y & y_y & z_y & 0 \\    x_z & y_z & z_z & 0...

Abstract This is a computer graphics and programming story. While ago, I wrote about a matrix that is generated by a OpenGL function, gluLookat. ( http://shitohichiumaya.blogspot.de/2011/01/what-matrix-glulookat-generates-1.html ) Why did I write that story? Because I was developing a renderer that is independent from the OpenGL, but some point, I need to overlay OpenGL generated image and our renderer's image. I think not a lot of people need to write a own renderer, therefore, this blog entry maybe not so much refereed. However, there are quite a lot of people visited this page. I don't know what kind of people visiting this pages, I am a bit surprised. This time, our customer has a special system and they need to compute the inverse of the OpenGL functions. Their input is OpenGL's Projection/Modelview matrix and the output is camera parameters. This is more special cases, but, the problem itself may be interested in. Introduction There are interesting visualization ...

Conclusion I think it is not exactly correct that the normal transform is an inverse transpose. A normal vector is not a column vector, the transform is a left inverse as Equation 5 shows. The transpose comes from the implementation of the vector class. This should be a left inverse as a rigorous formulation. I don't say the references [2,3,6] are incorrect, but I prefer the references [4,5] since they distinguish these two kind of vectors. However, the references [2,3,6] have an easier to understand explanations. Especially, [2,3,6] shows why normal transformation doesn't work with applying A intuitively. I like these explanations more. Here I would like to highlight the difference of a usual vector and a normal vector, they are a column vector or a row vector. By the way, contravariant vector means doesn't-change-vector, but, its representation can change. It seems I should explain the difference of these vectors are related with inner product. I would like to stud...

Why does it transpose of inverse matrix? In the last section, we know what is the transformation matrix for the normals. I will write it here again. But this is not a transpose of inverse matrix. It's just a left inverse matrix. Because a normal vector is a row vector, this is the correct notation. But, in a graphics library, i.e. OpenGL, we usually don't distinguish row vectors and column vectors in the computer memory. Moreover, we also don't distinguish points and vectors in the computer memory, they are usually length three (or four) of array. A usual vector is a column vector, and since we don't distinguish normal vectors and usual vectors, they are all treated as column vectors. But actually they are different. This point, you can see the difference. To make it normal vector, we need to transpose a normal vector to make it a column vector.  The transpose of Equation 5 is: Now you see why most of the books said normal transformation is a transpose of inve...

Let coordinate \Sigma 's the origin O , coordinate \Sigma' 's origin O' , then we can think about the coordinates of O'O in the coordinate system \Sigma' that is represented as: I would like to have a comment of this why this matters. In some graphics system, e.g., OpenGL, we can move or distort objects by applying transformation matrix.  This looks like the object is moved or is distorted. This interpretation is possible, but, here we did not move the objects, but we changed the coordinate system. This is rather how we interpret the result. I prefer this interpretation since we can think changing coordinate system without objects. You can still think the applying transformation matrix is an applying an operator, but, object is just one subject to operate. I would like to think rather about operator itself. Then we can concentrate operator itself, the transformation matrix itself. In Figure 3, there is a point P , this point actually doesn't move in ...

Normal vector as a covariant vector This third explanation includes how to transform normal vectors. This explanation is based on one of my favorite book[5] by Sugihara. This explanation is a bit formal and less intuitive compare to the first and the second explanation. If you are not interested in a formal explanation, you can skip this section. First I would like to introduce the affine transformation. Affine transformation is one of linear transformations. This is quite often used in computer graphics area. Affine transformation transform a line to a line and keep the ratio on a line. If we include a degenerated case, a triangle is always transformed into a triangle. Assume a representation of a three-dimensional affine transformation is a 3x3 matrix. This transformation is a transformation between two coodinate systems. Therefore, I think an object deformation by the transformation is a secondary effect. As a result, we can deform an object. However, this is a transformatio...

Normal defined by inner product Let's think about a normal is defined by inner product. As you see, the normal is actually in a row space instead of in a column space. A coordinate transformation matrix M is a transformation of a column vector. Therefore, the following equation is not defined. This is my second explanation what is the difference between normal vector and usual vector.

Normal vector as a perpendicular vector of the surface tangent vectors Normal vector has the same direction to the cross product of two tangent vectors of a surface.  Figure 2 shows the tangent vectors are correctly transformed by the matrix M that magnifies only x direction. However, their cross product is not necessary to the same as the transformation of normal by the matrix M . Figure 2. The normal vector n is a cross product of tangent vectors u and v . Tangent vectors are linear to M , but not for the normal vector. In short, tangent vectors u , v can be transformed by M , but their cross product is not. In general, Are you convinced this is the reason distinguishing a usual vector and a normal vector? If you think about the x component of the cross product,  u y v z - u z v y , this is not linear. Therefore, a linear transformation cannot transform this. This is my first explanation.

The difference between normal vectors and usual vectors Actually, normal vectors and usual vectors are different type of vectors. A usual vector itself defines direction or position. Let's recall how the normal vector is defined. A normal vector is not defined by itself. First a face is defined using usual vectors, then, the normal vector is defined according to this surface. The definition of the normal vector is a unit vector that perpendicular to the surface, in other words, it is defined as the inner product of surface tangent vector and the normal vector is 0. We should think these vectors are different vectors when they follow a transformation. As we saw in the Figure 1, when we transform a normal vector, the meaning of normal is lost. If we consider this vector is a usual vector, twice large x component is still fine, just it is not normal anymore. In this sense, normal vector has a special meaning, that should be always perpendicular to the surface. The following thr...

What is the problem? What is the transformation matrix of normal? We use transformation matrices every day when we move objects in a computer. Current the state of the art DCC (digital contents creation) software usually represents objects with triangles or polygons.  Each vertex of the triangles or polygons usually has its coordinates.  When we rotate or move each vertex, we apply a transformation matrix on each vertex. A vertex is usually three dimensional vector in computer graphics. We can define a normal vector for each triangle. A normal vector points out to which direction a triangle face is oriented. This normal vector is also a three dimensional vector. In a 3D computer graphics system, normal vectors are important since we need these normal vectors to compute how bright the surfaces are. Because a usual vector can be transformed by a matrix, it seems straightforward to use the same matrix to transform a normal vector. However, this fails. But why? The article i...

Abstract Several books[2,6] explained the normal vector transformation matrix is $(M^{-1})^{T}$ . I always forget this formula. This time I understand it a bit in three different ways, so I will write them down here. Introduction Assume matrix M is applied to a vector where the M is a coordinate transformation matrix. For example, M could be a translation, rotation, scaling, and so forth. To transform a position vector, we can just multiply this matrix M . However, we may fail when we transform a normal vector by just multiplying the matrix M . Several books mentioned normal transformation matrix should be $ (M^{-1})^{T}$ [2,6]. In this article, I would like to mention about the following three issues: What is the problem? What is the transformation matrix of normal? Why may multiplying matrix M fail? Why is it an inverse of transpose of matrix M ? (The references will be shown up at the end of this series.) Next time, I would like to talk about the problem.

gluLookAt implementation in python Here is my gluLookAt matrix computation code by python (numpy) and its test code. This code cares rather comprehensiveness than efficiency, therefore, there are matrix multiplication and matrix transpose. You can remove this if you want to have more efficient code. The test compares the result of this implementation and the OpenGL's computation result with less than 1.0e-6 error. As far as I tested, the result is in the error threshold in Ubuntu 10.04, python 2.6, glx version 1.4 ----- # Copyright (C) 2010-2011 H. Yamauchi # under New (3-clause) BSD license # # get lookat matrix (gluLookAt compatible matrix) python, numpy # # \param[in]  _eye    eye point # \param[in]  _lookat lookat point # \param[in]  _up     up vector # \return 4x4 gluLookAt matrix def getLookAtMatrix(_eye, _lookat, _up):   ez = _eye - _lookat   ez = ez / numpy.linalg.norm(ez)   ex = numpy.cross(_up,...

gluLookAt matrix (why Gimbal lock happens?)  The camera posture is represented by a rotation matrix, and camera position is represented by a translation matrix. Because in OpenGL, the default camera position and posture are defined and we move this camera around in the program. The camera is at origin, look into the Z axis minus direction, and up direction is Y axis plus direction. A rotation matrix and a translation matrix look like the following. where, x : camera basis X axis, ``right'' in the figure, ``side'' in the Mesa program in the last blog entry y : camera basis Y axis, ``up'' in the figure, ``up'' in the Mesa program in the last blog entry z : camera basis z axis, ``Z'' in the figure, ``-forward'' in the Mesa program in the last blog entry e : camera position, ``eye'' in the figure, ``eyex, eyey, eyez'' in the Mesa program in the last blog entry In the Mesa program, eye , lookat , up are given to t...

gluLookAt function I use OpenGL API to draw a 3D scene. This API is often used in games. If you only use OpenGL world, you usually don't need to know what kind of matrices are created underneath the API. However, if you want to overlay the OpenGL scene, for instance, if you write a own renderer and still want to render the scene with OpenGL and own, it is sometimes useful to know what kind of matrices are computed in the OpenGL API. Therefore, this blog entry is quite specific. If you are not interested in such story, see you next time. I would like to talk about a matrix that is computed by gluLookAt function. gluLookAt puts a camera with a position and camera posture, e.g., which is up direction. Recently, I need to know how to generate this matrix by my own, however, I could not find a good page with a good explanation. At the end, I looked into the source code of Mesa 7.5.1 (Mesa-7.5.1/src/glu/sgi/libutil/project.c), OpenGL software implementation. The code is quite r...

This issue of Communications of the ACM has an article by Igarashi, Takeo. http://doi.acm.org/10.1145/1785414.1785436 I liked this article since I can see the overview of his research history. If you allow my presumption to mention his research: his research has full of fun demos, is practical, is developing solid base technologies. They are great. I think people like his full of dream research. Yes, I also like that. But, I especially like the technology oriented papers from his team which I feel less high-lighted from the community. One example is SmoothTeddy. It seems that work is not treated well like as Teddy. But, I think SmoothTeddy has a significant technological contribution for mesh optimization for interactive sketching systems. I like his technical contribution like those papers.

While ago, I have a blog about a template error. I finally finished it today. I made a Sierpinski tetrahedron generator. http://en.wikipedia.org/wiki/Sierpinski_triangle This time, I would like to learn OpenMesh ( http://www.openmesh.org ). I have my own mesh library, but usually I would like to solve a problem, so if there is a well developed open source library, I think it is better to switch to it. Figure 1 shows how to get consistent list of vertices by OpenMesh. Here, "consistent" means, I always would like to have the vertices indexed as in the Figure 1, up. For example, if O_1 and O_3 are exchanged, I have a problem to create a Sierpinski tetrahedron. Figure 1: tetrahedron configuration Figure2: creation rule 1 Figure 3: creation rule2 Actually, I don't know what is the best way to do that, I depends on the OpenMesh's face halfedge circulator. If anybody knows better way, please give me a comment. In my method, first I need to guarantee the input is a te...

Page 50: Stratified I was not familiar with the word `Stratified.' A Monte-Carlo method usually uses a pseudo random number sequence. A random number sequence usually should have a property, `unexpectedness,' this means a same number sequence like 1, 1, 1, 1, 1, ..., 1 could be also a part of random number sequence. This issue is well explained in Knuth's book (D.E.Knuth, The Art of Computer Programming, Vol.2 Random numbers). The sampling pattern could be the left of Figure 1.  Figure 1. Stratified example. Left example of random 12 samples (clumping problem), Right example of 2x2 stratified random 12 samples In Figure 1 left, some part is not well sampled, or some part is oversampled -- the sampling pattern is biased. This is clumping, a bad luck in a sense. Usually it is solved when the number of samples is increased, but, it computationally costs, means it takes more time to compute. To solve this clumping problem without losing performance, the sampling region is ...

Page 11: Light implementation In page 11, there is a paragraph about light implementation. Ideally, the performance of the light transport algorithm should depend only on what the scene represents, rather than the details of how it is represented. For example, consider a scene illuminated by a square area light source. If this light source is replaced with a 10 by 10 grid of point sources, then the visual results will be nearly identical. However, the performance of many light transport algorithms will be much worse in the second case. Why 10x10 points is more expensive? As an amateur, I think if the area light is sampled by 20x20, then it seems 10x10 points is cheaper. But this is just an example. It seems the author want to say if you tell the renderer more straight forward way, the rendere can usually optimize the process. It is better to say that an area light is an area light, not an area light is many points. (Thanks to D.S.) Page 14: unbias and consistent This page mentio...

Here is my Japanese translation, part 3/3. Robust Monte Carlo Methods for Light Transport Simulation 光輸送シミュレーションのためのロバストなモンテカルロ法 著者(Author): Eric Veach Translated by Lx=d HY 1.1.2 輸送モデルに関する仮定 光輸送アルゴリズムは全ての光の振舞いをシミュレートするのではない．なぜなら，ほとんどの応用でそれは必要でないからである．(脚注: 厳密に言えばそれは可能ですらない．なぜなら全ての物理学が知られているわけではないからである．しかし，光と物体の相互作用は物理学の中でも最もよく研究され，発展しているものの一つである．そして実質的に全ての観測される現象は，かなり高い精度で計算，予測できる．[Feynman 1985] したがって，コンピュータグラフィックスの目的としては我々はこれらの物理法則は，ほぼ完全に理解されていると仮定してよい．) グラフィックスの立場に立つと，物理学における光学の知識を使うことは，写実的な画像を生成する手法として，最も適していると考えられる．どのようなシーンを生成するかに対して，我々はどの光学的効果が重要なのかを判断し，そしてそれをシミュレートできるアルゴリズムを選ぶことになる． 我々の論文においては幾何光学モデルを基礎とする．このモデルでは，光は物体の表面からのみ放射，拡散，吸収される．そして，これらの面間では光は直進すると仮定する．したがって，雲や煙のような participate 物体(吸収・拡散を行う物体)や，あるいは屈折率が連続的に変化する物体(たとえば，熱せられた気体)は考慮しない．我々はまた，波長や量子力学的効果に依存する(たとえば回折や蛍光)光学効果のほとんどを無視する．特に，我々は光線間の相互作用の可能性を排除する．すなわち，光は完全に非干渉であることを仮定する． 通常のレンダリングでは，ここで我々が無視した効果はそれほど重要ではない．幾何光学は我々が日常で目にするものをほとんど全て，適切に高い精度でモデル化する．このため，実質的には，ほとんどのグラフィックスでの光輸送アルゴリズムはここで行なった仮定と同様の仮定を行っている．こ...