Here is the second assignment @ Link.
The output looks like this:
I lazily kept postponing working on this assignment and somewhat slogged in the last few days. The professor kept postponing the last date for submitting the assignment so that helped. My purpose of posting the solutions is not to help future students but with an intension of forcing changes in assignments and updation of questions.
I don't know if the gimbal lock problem is there in this solution or not. I thought using quaternions solves the problem. If you think there is a problem still plz contact me.
I didn't solve this assignment completely. Anyone with basic programming skills can modify the code and make it work according to the expectations of the professor.
Here it is:
/*
* ellipsoid.c
*
* Created on: Oct 27, 2013
* Author: kamath
*/
#include<GL/glut.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<unistd.h>
#include "quater.h"
#define Pi 3.141
static GLdouble viewer[]={80,80,50};
float phi=0,timestamp,alpha1,beta1,gamma1;
FILE *fp;
void DrawEllipsoid(unsigned int uiStacks, unsigned int uiSlices, float fA, float fB, float fC)
{
float s,t,tStep = (Pi) / (float)uiSlices;
float sStep = (Pi) / (float)uiStacks;
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
for(t = -Pi/2; t <= (Pi/2)+.0001; t += tStep)
{
glBegin(GL_TRIANGLE_STRIP);
for(s = -Pi; s <= Pi+.0001; s += sStep)
{
glVertex3f(fA * cos(t) * cos(s), fB * cos(t) * sin(s), fC * sin(t));
glVertex3f(fA * cos(t+tStep) * cos(s), fB * cos(t+tStep) * sin(s), fC * sin(t+tStep));
}
glEnd();
}
}
void display()
{float m2[4][4],m3[4][4],m4[4][4],m21[16],m31[16],m41[16];
float a1[3]={0,0,1.0},a2[3]={0,1.0,0},q[4],k[4],j[4];
int i;
alpha1=alpha1*3.141/180; //convert degrees into radians
beta1=beta1*3.141/180;
gamma1=gamma1*3.141/180;
axis_to_quat(a1,alpha1,q); //convert euler angles into quaternions
axis_to_quat(a2,beta1,k);
axis_to_quat(a1,gamma1,j);
//add_quats(q,k,dest)
build_rotmatrix(m2,q); //convert quaternions into matrices
build_rotmatrix(m3,k);
build_rotmatrix(m4,j);
for(i=0;i<16;i++) //convert matrix into column major order
{
m21[i]=m2[i%4][(int)floor(i/4)];
}
for(i=0;i<16;i++)
{
m31[i]=m3[i%4][(int)floor(i/4)];
}
for(i=0;i<16;i++)
{
m41[i]=m4[i%4][(int)floor(i/4)];
}
usleep(1000);
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
glColor3f(1,0,0);
glLoadIdentity();
gluLookAt(viewer[0],viewer[1],viewer[2],0,0,0,0,1,0);
glLoadIdentity();
glPushMatrix();
glRotatef(45,1,1,1);
glBegin(GL_LINES);
glVertex3f(0,0,-10);
glVertex3f(0,0,10);
glVertex3f(10,0,0);
glVertex3f(-10,0,0);
glVertex3f(0,10,0);
glVertex3f(0,-10,0);
glEnd();
glColor3f(0.8,0.8,0.4);
glTranslatef(30,30,0);
glMultMatrixf(m41);//rotate wrt Z
glMultMatrixf(m31);//rotate wrt x
glMultMatrixf(m21);//rotate wrt z
DrawEllipsoid(10, 10, 12,6,24);
glFlush();
glPopMatrix();
glutSwapBuffers();
}
void init()
{
glClearColor(1.0,1.0,1.0,0.0);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(-80,80,-80,80,-80,80);
glMatrixMode(GL_MODELVIEW);
}
void idle()
{// read file and get angles
if(fscanf (fp,"%f %f %f %f\n",×tamp,&alpha1,&beta1,&gamma1)==EOF)
{
fclose(fp);
printf("end of input");
exit(0);
}
glutPostRedisplay();
}
int main(int argc,char *argv[])
{
fp = fopen ("el1_sh.txt", "r");
if (fp == NULL)
{
printf ("Error opening file");
exit(1);
}
glutInit(&argc,argv);
glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB);
glutInitWindowSize(500,500);
glutInitWindowPosition(200,0);
glutCreateWindow("ellipsoid");
glutDisplayFunc(display);
init();
glEnable(GL_DEPTH_TEST);
glutIdleFunc(idle);
glutMainLoop();
return 0;
}
This is a .h file you need to include: (Contains code to work with quaternions)
#include <math.h>
//#include "trackball.h"
/*
* This size should really be based on the distance from the center of
* rotation to the point on the object underneath the mouse. That
* point would then track the mouse as closely as possible. This is a
* simple example, though, so that is left as an Exercise for the
* Programmer.
*/
#define TRACKBALLSIZE (0.8)
/*
* Local function prototypes (not defined in trackball.h)
*/
static float tb_project_to_sphere(float, float, float);
static void normalize_quat(float [4]);
void
vzero(float *v)
{
v[0] = 0.0;
v[1] = 0.0;
v[2] = 0.0;
}
void
vset(float *v, float x, float y, float z)
{
v[0] = x;
v[1] = y;
v[2] = z;
}
void
vsub(const float *src1, const float *src2, float *dst)
{
dst[0] = src1[0] - src2[0];
dst[1] = src1[1] - src2[1];
dst[2] = src1[2] - src2[2];
}
void
vcopy(const float *v1, float *v2)
{
register int i;
for (i = 0 ; i < 3 ; i++)
v2[i] = v1[i];
}
void
vcross(const float *v1, const float *v2, float *cross)
{
float temp[3];
temp[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
temp[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
temp[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
vcopy(temp, cross);
}
float
vlength(const float *v)
{
return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}
void
vscale(float *v, float div)
{
v[0] *= div;
v[1] *= div;
v[2] *= div;
}
void
vnormal(float *v)
{
vscale(v,1.0/vlength(v));
}
float
vdot(const float *v1, const float *v2)
{
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
void
vadd(const float *src1, const float *src2, float *dst)
{
dst[0] = src1[0] + src2[0];
dst[1] = src1[1] + src2[1];
dst[2] = src1[2] + src2[2];
}
/*
* Given an axis and angle, compute quaternion.
*/
void
axis_to_quat(float a[3], float phi, float q[4])
{
vnormal(a);
vcopy(a,q);
vscale(q,sin(phi/2.0));
q[3] = cos(phi/2.0);
}
/*
* Ok, simulate a track-ball. Project the points onto the virtual
* trackball, then figure out the axis of rotation, which is the cross
* product of P1 P2 and O P1 (O is the center of the ball, 0,0,0)
* Note: This is a deformed trackball-- is a trackball in the center,
* but is deformed into a hyperbolic sheet of rotation away from the
* center. This particular function was chosen after trying out
* several variations.
*
* It is assumed that the arguments to this routine are in the range
* (-1.0 ... 1.0)
*/
void
trackball(float q[4], float p1x, float p1y, float p2x, float p2y)
{
float a[3]; /* Axis of rotation */
float phi; /* how much to rotate about axis */
float p1[3], p2[3], d[3];
float t;
if (p1x == p2x && p1y == p2y) {
/* Zero rotation */
vzero(q);
q[3] = 1.0;
return;
}
/*
* First, figure out z-coordinates for projection of P1 and P2 to
* deformed sphere
*/
vset(p1,p1x,p1y,tb_project_to_sphere(TRACKBALLSIZE,p1x,p1y));
vset(p2,p2x,p2y,tb_project_to_sphere(TRACKBALLSIZE,p2x,p2y));
/*
* Now, we want the cross product of P1 and P2
*/
vcross(p2,p1,a);
/*
* Figure out how much to rotate around that axis.
*/
vsub(p1,p2,d);
t = vlength(d) / (2.0*TRACKBALLSIZE);
/*
* Avoid problems with out-of-control values...
*/
if (t > 1.0) t = 1.0;
if (t < -1.0) t = -1.0;
phi = 2.0 * asin(t);
axis_to_quat(a,phi,q);
}
/*
* Project an x,y pair onto a sphere of radius r OR a hyperbolic sheet
* if we are away from the center of the sphere.
*/
static float
tb_project_to_sphere(float r, float x, float y)
{
float d, t, z;
d = sqrt(x*x + y*y);
if (d < r * 0.70710678118654752440) { /* Inside sphere */
z = sqrt(r*r - d*d);
} else { /* On hyperbola */
t = r / 1.41421356237309504880;
z = t*t / d;
}
return z;
}
/*
* Given two rotations, e1 and e2, expressed as quaternion rotations,
* figure out the equivalent single rotation and stuff it into dest.
*
* This routine also normalizes the result every RENORMCOUNT times it is
* called, to keep error from creeping in.
*
* NOTE: This routine is written so that q1 or q2 may be the same
* as dest (or each other).
*/
#define RENORMCOUNT 97
void
add_quats(float q1[4], float q2[4], float dest[4])
{
static int count=0;
float t1[4], t2[4], t3[4];
float tf[4];
vcopy(q1,t1);
vscale(t1,q2[3]);
vcopy(q2,t2);
vscale(t2,q1[3]);
vcross(q2,q1,t3);
vadd(t1,t2,tf);
vadd(t3,tf,tf);
tf[3] = q1[3] * q2[3] - vdot(q1,q2);
dest[0] = tf[0];
dest[1] = tf[1];
dest[2] = tf[2];
dest[3] = tf[3];
if (++count > RENORMCOUNT) {
count = 0;
normalize_quat(dest);
}
}
/*
* Quaternions always obey: a^2 + b^2 + c^2 + d^2 = 1.0
* If they don't add up to 1.0, dividing by their magnitued will
* renormalize them.
*
* Note: See the following for more information on quaternions:
*
* - Shoemake, K., Animating rotation with quaternion curves, Computer
* Graphics 19, No 3 (Proc. SIGGRAPH'85), 245-254, 1985.
* - Pletinckx, D., Quaternion calculus as a basic tool in computer
* graphics, The Visual Computer 5, 2-13, 1989.
*/
static void
normalize_quat(float q[4])
{
int i;
float mag;
mag = (q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
for (i = 0; i < 4; i++) q[i] /= mag;
}
/*
* Build a rotation matrix, given a quaternion rotation.
*
*/
void
build_rotmatrix(float m[4][4], float q[4])
{
m[0][0] = 1.0 - 2.0 * (q[1] * q[1] + q[2] * q[2]);
m[0][1] = 2.0 * (q[0] * q[1] - q[2] * q[3]);
m[0][2] = 2.0 * (q[2] * q[0] + q[1] * q[3]);
m[0][3] = 0.0;
m[1][0] = 2.0 * (q[0] * q[1] + q[2] * q[3]);
m[1][1]= 1.0 - 2.0 * (q[2] * q[2] + q[0] * q[0]);
m[1][2] = 2.0 * (q[1] * q[2] - q[0] * q[3]);
m[1][3] = 0.0;
m[2][0] = 2.0 * (q[2] * q[0] - q[1] * q[3]);
m[2][1] = 2.0 * (q[1] * q[2] + q[0] * q[3]);
m[2][2] = 1.0 - 2.0 * (q[1] * q[1] + q[0] * q[0]);
m[2][3] = 0.0;
m[3][0] = 0.0;
m[3][1] = 0.0;
m[3][2] = 0.0;
m[3][3] = 1.0;
}
This is a text file el1_sh.txt that is the input to the rotation functions, so u need to put it inside the project folder. I will soon update the post to include a few more comments to explain the program.
The output looks like this:
I lazily kept postponing working on this assignment and somewhat slogged in the last few days. The professor kept postponing the last date for submitting the assignment so that helped. My purpose of posting the solutions is not to help future students but with an intension of forcing changes in assignments and updation of questions.
I don't know if the gimbal lock problem is there in this solution or not. I thought using quaternions solves the problem. If you think there is a problem still plz contact me.
I didn't solve this assignment completely. Anyone with basic programming skills can modify the code and make it work according to the expectations of the professor.
Here it is:
/*
* ellipsoid.c
*
* Created on: Oct 27, 2013
* Author: kamath
*/
#include<GL/glut.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<unistd.h>
#include "quater.h"
#define Pi 3.141
static GLdouble viewer[]={80,80,50};
float phi=0,timestamp,alpha1,beta1,gamma1;
FILE *fp;
void DrawEllipsoid(unsigned int uiStacks, unsigned int uiSlices, float fA, float fB, float fC)
{
float s,t,tStep = (Pi) / (float)uiSlices;
float sStep = (Pi) / (float)uiStacks;
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
for(t = -Pi/2; t <= (Pi/2)+.0001; t += tStep)
{
glBegin(GL_TRIANGLE_STRIP);
for(s = -Pi; s <= Pi+.0001; s += sStep)
{
glVertex3f(fA * cos(t) * cos(s), fB * cos(t) * sin(s), fC * sin(t));
glVertex3f(fA * cos(t+tStep) * cos(s), fB * cos(t+tStep) * sin(s), fC * sin(t+tStep));
}
glEnd();
}
}
void display()
{float m2[4][4],m3[4][4],m4[4][4],m21[16],m31[16],m41[16];
float a1[3]={0,0,1.0},a2[3]={0,1.0,0},q[4],k[4],j[4];
int i;
alpha1=alpha1*3.141/180; //convert degrees into radians
beta1=beta1*3.141/180;
gamma1=gamma1*3.141/180;
axis_to_quat(a1,alpha1,q); //convert euler angles into quaternions
axis_to_quat(a2,beta1,k);
axis_to_quat(a1,gamma1,j);
//add_quats(q,k,dest)
build_rotmatrix(m2,q); //convert quaternions into matrices
build_rotmatrix(m3,k);
build_rotmatrix(m4,j);
for(i=0;i<16;i++) //convert matrix into column major order
{
m21[i]=m2[i%4][(int)floor(i/4)];
}
for(i=0;i<16;i++)
{
m31[i]=m3[i%4][(int)floor(i/4)];
}
for(i=0;i<16;i++)
{
m41[i]=m4[i%4][(int)floor(i/4)];
}
usleep(1000);
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
glColor3f(1,0,0);
glLoadIdentity();
gluLookAt(viewer[0],viewer[1],viewer[2],0,0,0,0,1,0);
glLoadIdentity();
glPushMatrix();
glRotatef(45,1,1,1);
glBegin(GL_LINES);
glVertex3f(0,0,-10);
glVertex3f(0,0,10);
glVertex3f(10,0,0);
glVertex3f(-10,0,0);
glVertex3f(0,10,0);
glVertex3f(0,-10,0);
glEnd();
glColor3f(0.8,0.8,0.4);
glTranslatef(30,30,0);
glMultMatrixf(m41);//rotate wrt Z
glMultMatrixf(m31);//rotate wrt x
glMultMatrixf(m21);//rotate wrt z
DrawEllipsoid(10, 10, 12,6,24);
glFlush();
glPopMatrix();
glutSwapBuffers();
}
void init()
{
glClearColor(1.0,1.0,1.0,0.0);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glOrtho(-80,80,-80,80,-80,80);
glMatrixMode(GL_MODELVIEW);
}
void idle()
{// read file and get angles
if(fscanf (fp,"%f %f %f %f\n",×tamp,&alpha1,&beta1,&gamma1)==EOF)
{
fclose(fp);
printf("end of input");
exit(0);
}
glutPostRedisplay();
}
int main(int argc,char *argv[])
{
fp = fopen ("el1_sh.txt", "r");
if (fp == NULL)
{
printf ("Error opening file");
exit(1);
}
glutInit(&argc,argv);
glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB);
glutInitWindowSize(500,500);
glutInitWindowPosition(200,0);
glutCreateWindow("ellipsoid");
glutDisplayFunc(display);
init();
glEnable(GL_DEPTH_TEST);
glutIdleFunc(idle);
glutMainLoop();
return 0;
}
This is a .h file you need to include: (Contains code to work with quaternions)
#include <math.h>
//#include "trackball.h"
/*
* This size should really be based on the distance from the center of
* rotation to the point on the object underneath the mouse. That
* point would then track the mouse as closely as possible. This is a
* simple example, though, so that is left as an Exercise for the
* Programmer.
*/
#define TRACKBALLSIZE (0.8)
/*
* Local function prototypes (not defined in trackball.h)
*/
static float tb_project_to_sphere(float, float, float);
static void normalize_quat(float [4]);
void
vzero(float *v)
{
v[0] = 0.0;
v[1] = 0.0;
v[2] = 0.0;
}
void
vset(float *v, float x, float y, float z)
{
v[0] = x;
v[1] = y;
v[2] = z;
}
void
vsub(const float *src1, const float *src2, float *dst)
{
dst[0] = src1[0] - src2[0];
dst[1] = src1[1] - src2[1];
dst[2] = src1[2] - src2[2];
}
void
vcopy(const float *v1, float *v2)
{
register int i;
for (i = 0 ; i < 3 ; i++)
v2[i] = v1[i];
}
void
vcross(const float *v1, const float *v2, float *cross)
{
float temp[3];
temp[0] = (v1[1] * v2[2]) - (v1[2] * v2[1]);
temp[1] = (v1[2] * v2[0]) - (v1[0] * v2[2]);
temp[2] = (v1[0] * v2[1]) - (v1[1] * v2[0]);
vcopy(temp, cross);
}
float
vlength(const float *v)
{
return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}
void
vscale(float *v, float div)
{
v[0] *= div;
v[1] *= div;
v[2] *= div;
}
void
vnormal(float *v)
{
vscale(v,1.0/vlength(v));
}
float
vdot(const float *v1, const float *v2)
{
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2];
}
void
vadd(const float *src1, const float *src2, float *dst)
{
dst[0] = src1[0] + src2[0];
dst[1] = src1[1] + src2[1];
dst[2] = src1[2] + src2[2];
}
/*
* Given an axis and angle, compute quaternion.
*/
void
axis_to_quat(float a[3], float phi, float q[4])
{
vnormal(a);
vcopy(a,q);
vscale(q,sin(phi/2.0));
q[3] = cos(phi/2.0);
}
/*
* Ok, simulate a track-ball. Project the points onto the virtual
* trackball, then figure out the axis of rotation, which is the cross
* product of P1 P2 and O P1 (O is the center of the ball, 0,0,0)
* Note: This is a deformed trackball-- is a trackball in the center,
* but is deformed into a hyperbolic sheet of rotation away from the
* center. This particular function was chosen after trying out
* several variations.
*
* It is assumed that the arguments to this routine are in the range
* (-1.0 ... 1.0)
*/
void
trackball(float q[4], float p1x, float p1y, float p2x, float p2y)
{
float a[3]; /* Axis of rotation */
float phi; /* how much to rotate about axis */
float p1[3], p2[3], d[3];
float t;
if (p1x == p2x && p1y == p2y) {
/* Zero rotation */
vzero(q);
q[3] = 1.0;
return;
}
/*
* First, figure out z-coordinates for projection of P1 and P2 to
* deformed sphere
*/
vset(p1,p1x,p1y,tb_project_to_sphere(TRACKBALLSIZE,p1x,p1y));
vset(p2,p2x,p2y,tb_project_to_sphere(TRACKBALLSIZE,p2x,p2y));
/*
* Now, we want the cross product of P1 and P2
*/
vcross(p2,p1,a);
/*
* Figure out how much to rotate around that axis.
*/
vsub(p1,p2,d);
t = vlength(d) / (2.0*TRACKBALLSIZE);
/*
* Avoid problems with out-of-control values...
*/
if (t > 1.0) t = 1.0;
if (t < -1.0) t = -1.0;
phi = 2.0 * asin(t);
axis_to_quat(a,phi,q);
}
/*
* Project an x,y pair onto a sphere of radius r OR a hyperbolic sheet
* if we are away from the center of the sphere.
*/
static float
tb_project_to_sphere(float r, float x, float y)
{
float d, t, z;
d = sqrt(x*x + y*y);
if (d < r * 0.70710678118654752440) { /* Inside sphere */
z = sqrt(r*r - d*d);
} else { /* On hyperbola */
t = r / 1.41421356237309504880;
z = t*t / d;
}
return z;
}
/*
* Given two rotations, e1 and e2, expressed as quaternion rotations,
* figure out the equivalent single rotation and stuff it into dest.
*
* This routine also normalizes the result every RENORMCOUNT times it is
* called, to keep error from creeping in.
*
* NOTE: This routine is written so that q1 or q2 may be the same
* as dest (or each other).
*/
#define RENORMCOUNT 97
void
add_quats(float q1[4], float q2[4], float dest[4])
{
static int count=0;
float t1[4], t2[4], t3[4];
float tf[4];
vcopy(q1,t1);
vscale(t1,q2[3]);
vcopy(q2,t2);
vscale(t2,q1[3]);
vcross(q2,q1,t3);
vadd(t1,t2,tf);
vadd(t3,tf,tf);
tf[3] = q1[3] * q2[3] - vdot(q1,q2);
dest[0] = tf[0];
dest[1] = tf[1];
dest[2] = tf[2];
dest[3] = tf[3];
if (++count > RENORMCOUNT) {
count = 0;
normalize_quat(dest);
}
}
/*
* Quaternions always obey: a^2 + b^2 + c^2 + d^2 = 1.0
* If they don't add up to 1.0, dividing by their magnitued will
* renormalize them.
*
* Note: See the following for more information on quaternions:
*
* - Shoemake, K., Animating rotation with quaternion curves, Computer
* Graphics 19, No 3 (Proc. SIGGRAPH'85), 245-254, 1985.
* - Pletinckx, D., Quaternion calculus as a basic tool in computer
* graphics, The Visual Computer 5, 2-13, 1989.
*/
static void
normalize_quat(float q[4])
{
int i;
float mag;
mag = (q[0]*q[0] + q[1]*q[1] + q[2]*q[2] + q[3]*q[3]);
for (i = 0; i < 4; i++) q[i] /= mag;
}
/*
* Build a rotation matrix, given a quaternion rotation.
*
*/
void
build_rotmatrix(float m[4][4], float q[4])
{
m[0][0] = 1.0 - 2.0 * (q[1] * q[1] + q[2] * q[2]);
m[0][1] = 2.0 * (q[0] * q[1] - q[2] * q[3]);
m[0][2] = 2.0 * (q[2] * q[0] + q[1] * q[3]);
m[0][3] = 0.0;
m[1][0] = 2.0 * (q[0] * q[1] + q[2] * q[3]);
m[1][1]= 1.0 - 2.0 * (q[2] * q[2] + q[0] * q[0]);
m[1][2] = 2.0 * (q[1] * q[2] - q[0] * q[3]);
m[1][3] = 0.0;
m[2][0] = 2.0 * (q[2] * q[0] - q[1] * q[3]);
m[2][1] = 2.0 * (q[1] * q[2] + q[0] * q[3]);
m[2][2] = 1.0 - 2.0 * (q[1] * q[1] + q[0] * q[0]);
m[2][3] = 0.0;
m[3][0] = 0.0;
m[3][1] = 0.0;
m[3][2] = 0.0;
m[3][3] = 1.0;
}
This is a text file el1_sh.txt that is the input to the rotation functions, so u need to put it inside the project folder. I will soon update the post to include a few more comments to explain the program.
No comments:
Post a Comment