Calculating normals of closed objects

Hi, I have a problem with the calculation of the normals of an object. I followed the terrain tutorial from VideoTutorialsRock (see http://www.videotutorialsrock.com/opengl_tutorial/terrain/home.php)) and I wanted to close the terrain surface. In this I have succeeded, however the shading of the surface does not seem correct. At the point where the edges of the original surface are brought together you clearly see the edge.
I think this problem is due to the calculation of the normals, could anyone help?

Thanks!!

(You can see a picture of what I mean on http://users.telenet.be/lalaland/wouter/terrein.jpg )

#include <iostream>
#include <stdlib.h>
#include <math.h>

#ifdef __APPLE__
#include <OpenGL/OpenGL.h>
#include <GLUT/glut.h>
#else
#include <GL/glut.h>
#endif

#include "imageloader.h"
#include "vec3f.h"

#define PI 3.14159265

using namespace std;

//Represents a terrain, by storing a set of heights and normals at 2D locations
class Terrain {
	private:
		int w; //Width
		int l; //Length
		float r; //Radius
		float*** pos; // xPos, yPos, zPos
		Vec3f** normals;
		bool extra;
		bool computedNormals; //Whether normals is up-to-date
	public:
		Terrain(int w2, int l2, int r2) {
			w = w2;
			l = l2;
			r = r2;

			pos = new float**[l];
			normals = new Vec3f*[l];
			for(int i = 0; i < l; i++) {
				pos[i] = new float*[w];
				normals[i] = new Vec3f[w];
				for (int j=0; j<w; j++)
					pos[i][j] = new float[3];
			}
			
			extra = false;
			computedNormals = false;

			for (int i = 0; i < l; i++) {
				for (int j = 0; j < w; j++) {
					pos[i][j][0] = r*sin(j*2*PI/(w-1));
					pos[i][j][1] = r*cos(j*2*PI/(w-1));
					pos[i][j][2] = i;
				}
			}
		}
		
		~Terrain() {
			for(int i = 0; i < l; i++) {
				for(int j = 0; j < w; j++) {
					delete[] pos[i][j];
				}
				delete[] pos[i];
				delete[] normals[i];
			}
			delete[] pos;
			delete[] normals;
		}
		
		int width() {
			return w;
		}
		
		int length() {
			return l;
		}

		int radius() {
			return r;
		}
		
		//Sets the height at (x, z) to h
		void setHeight(int x, int z, float h) {
			pos[z][x][0] = (r+h)*sin(x*2*PI/(w-1));
			pos[z][x][1] = (r+h)*cos(x*2*PI/(w-1));
			computedNormals = false;
		}
		
		float getXPos(int x, int z) {
			return pos[z][x][0];
		}

		float getYPos(int x, int z) {
			return pos[z][x][1];
		}

		float getZPos(int x, int z) {
			return pos[z][x][2];
		}

		//Returns the height at (x, z)
		float getHeight(int x, int z) {
			return sqrt(pow(pos[z][x][0],2) + pow(pos[z][x][1],2)) - r;
		}

		void switchExtra() {
			computedNormals = false;
			extra = !extra;
		}
		
		//Computes the normals, if they haven't been computed yet
		void computeNormals() {
			if (computedNormals) {
				return;
			}
			
			//Compute the rough version of the normals
			Vec3f** normals2 = new Vec3f*[l];
			for(int i = 0; i < l; i++) {
				normals2[i] = new Vec3f[w];
			}
			
			for(int z = 0; z < l; z++) {
				for(int x = 0; x < w; x++) {
					Vec3f sum(0.0f, 0.0f, 0.0f);
					
					Vec3f out;
					if (z > 0) {
						out = Vec3f(pos[z - 1][x][0] - pos[z][x][0], pos[z - 1][x][1] - pos[z][x][1], pos[z - 1][x][2] - pos[z][x][2]);
					}
					Vec3f in;
					if (z < l - 1) {
						in = Vec3f(pos[z + 1][x][0] - pos[z][x][0], pos[z + 1][x][1] - pos[z][x][1], pos[z + 1][x][2] - pos[z][x][2]);
					}
					Vec3f left;
					left = Vec3f(pos[z][(x - 1+w)%w][0] - pos[z][x][0], pos[z][(x - 1+w)%w][1] - pos[z][x][1], pos[z][(x - 1+w)%w][2] - pos[z][x][2]);
					Vec3f right;
					right = Vec3f(pos[z][(x + 1)%w][0] - pos[z][x][0], pos[z][(x + 1)%w][1] - pos[z][x][1], pos[z][(x + 1)%w][2] - pos[z][x][2]);
					
					if (z > 0) {
						sum += out.cross(left).normalize();
						sum += right.cross(out).normalize();
					}
					if (z < l - 1) {
						sum += left.cross(in).normalize();
						sum += in.cross(right).normalize();
					}
					
					normals2[z][x] = sum;
				}
			}
			
			//Smooth out the normals
			const float FALLOUT_RATIO = 0.5f;
			for(int z = 0; z < l; z++) {
				for(int x = 0; x < w; x++) {
					Vec3f sum = normals2[z][x];
					
						sum += normals2[z][(x - 1)%w] * FALLOUT_RATIO;
						sum += normals2[z][(x + 1)%w] * FALLOUT_RATIO;
					if (z > 0) {
						sum += normals2[z - 1][x] * FALLOUT_RATIO;
					}
					if (z < l - 1) {
						sum += normals2[z + 1][x] * FALLOUT_RATIO;
					}
					if(extra) {
						if (z > 0) {
							sum += normals2[z - 1][(x - 1+w)%w] * FALLOUT_RATIO;
							sum += normals2[z - 1][(x + 1)%w] * FALLOUT_RATIO;
						}
						if (z < l - 1) {
							sum += normals2[z + 1][(x - 1+w)%w] * FALLOUT_RATIO;
							sum += normals2[z + 1][(x + 1)%w] * FALLOUT_RATIO;
						}
					}
					
					if (sum.magnitude() == 0) {
						sum = Vec3f(0.0f, 1.0f, 0.0f);
					}
					normals[z][x] = sum;
				}
			}
			
			for(int i = 0; i < l; i++) {
				delete[] normals2[i];
			}
			delete[] normals2;
			
			computedNormals = true;
		}
		
		//Returns the normal at (x, z)
		Vec3f getNormal(int x, int z) {
			if (!computedNormals) {
				computeNormals();
			}
			return normals[z][x];
		}
};

//Loads a terrain from a heightmap.  The heights of the terrain range from
//-height / 2 to height / 2.
Terrain* loadTerrain(const char* filename1, float height, float rad) {
	Image* image1 = loadBMP(filename1);
	Terrain* t1 = new Terrain(image1->width, image1->height, rad);
	for(int y = 0; y < image1->height; y++) {
		for(int x = 0; x < image1->width; x++) {
			unsigned char color1 =
				(unsigned char)image1->pixels[3 * (y * image1->width + x)];
			float h = height * color1 / 255.0f;
			t1->setHeight(x, y, h);
		}
	}
	
	delete image1;
	t1->computeNormals();
	return t1;
}

float _angle = 0.0f;
float _angleInc = 0.1f;
float _radius = 40.0f;
float _height = 20.0f;
Terrain* _terrain;

void cleanup() {
	delete _terrain;
}

void handleKeypress(unsigned char key, int x, int y) {
	switch (key) {
		case 27: //Escape key
			cleanup();
			exit(0);
			break;
		case 97: // 'a'
			_terrain->switchExtra();
			_terrain->computeNormals();
			glutPostRedisplay();
			break;
		case 113: // 'q'
			_angleInc -= 0.05f;
			break;
		case 115: // 's'
			_angleInc += 0.05f;
			break;
	}
}

void initRendering() {
	glEnable(GL_DEPTH_TEST);
	glEnable(GL_COLOR_MATERIAL);
	glEnable(GL_LIGHTING);
	glEnable(GL_LIGHT0);
	glEnable(GL_NORMALIZE);
	glShadeModel(GL_SMOOTH);
}

void handleResize(int w, int h) {
	glViewport(0, 0, w, h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(45.0, (double)w / (double)h, 1.0, 200.0);
}

void drawScene() {
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

	glRotatef(90.0f, 0.0f, 1.0f, 0.0f);
	glRotatef(_angle, 0.0f, 0.0f, -1.0f);
	glTranslatef((_radius+_height)*cos(_angle*2*PI/360), (_radius+_height)*sin(_angle*2*PI/360), 0.0f);

	GLfloat ambientColor[] = {0.4f, 0.4f, 0.4f, 1.0f};
	glLightModelfv(GL_LIGHT_MODEL_AMBIENT, ambientColor);
	
	GLfloat lightColor0[] = {0.6f, 0.6f, 0.6f, 1.0f};
	GLfloat lightPos0[] = {-0.5f, 0.8f, 0.1f, 0.0f};
	glLightfv(GL_LIGHT0, GL_DIFFUSE, lightColor0);
	glLightfv(GL_LIGHT0, GL_POSITION, lightPos0);
	
	float scale = 2*_radius / max(_terrain->width() - 1, _terrain->length() - 1);
	glScalef(scale, scale, scale);
	glTranslatef(0.0f,
				 0.0f,
				 -(float)(_terrain->length() - 1) / 2);
	
	glColor3f(0.9f, 0.7f, 0.0f);
	for(int z = 0; z < _terrain->length() - 1; z++) {
		//Makes OpenGL draw a triangle at every three consecutive vertices
		glBegin(GL_TRIANGLE_STRIP);
		for(int x = 0; x < _terrain->width(); x++) {
//			float theta = x*2*PI/(_terrain->width()-1);

//			float amplitude = _radius+_terrain->getHeight(x, z);
			Vec3f normal = _terrain->getNormal(x, z);
			glNormal3f(normal[0], normal[1], normal[2]);
//			glVertex3f(amplitude*sin(theta), amplitude*cos(theta), z);
			glVertex3f(_terrain->getXPos(x,z),_terrain->getYPos(x,z),_terrain->getZPos(x,z));

//			amplitude = _radius+_terrain->getHeight(x, z+1);
			normal = _terrain->getNormal(x, z + 1);
			glNormal3f(normal[0], normal[1], normal[2]);
//			glVertex3f(amplitude*sin(theta), amplitude*cos(theta), z+1);
			glVertex3f(_terrain->getXPos(x,z+1),_terrain->getYPos(x,z+1),_terrain->getZPos(x,z+1));
		}
		glEnd();
	}
	
	glutSwapBuffers();
}

void update(int value) {
	_angle += _angleInc;
	if (_angle > 360) {
		_angle -= 360;
	}
	
	glutPostRedisplay();
	glutTimerFunc(25, update, 0);
}

int main(int argc, char** argv) {
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glutInitWindowSize(400, 400);
	
	glutCreateWindow("Terrain - videotutorialsrock.com");
	initRendering();
	
	_terrain = loadTerrain("test5.bmp", _height, _radius);
	
	glutDisplayFunc(drawScene);
	glutKeyboardFunc(handleKeypress);
	glutReshapeFunc(handleResize);
	glutTimerFunc(25, update, 0);
	
	glutMainLoop();
	return 0;
}

For the additional headers, see http://www.videotutorialsrock.com/opengl_tutorial/terrain/home.php

Not sure but you might not be calculating the normals right
Here’s what I do

void CSurface::UpdateNormals()
{
	
	map<CVector3,CVector3> nmap;
	map<CVector3,CVector3> tmap;
	map<CVector3,CVector3> bmap;
	for(int i=0;i<this->tris_count;i++)
	{
		UINT v1 = tris[i*3+0];
		UINT v2 = tris[i*3+1];
		UINT v3 = tris[i*3+2];

		CVector3 vec1 = CVector3(vert_coords[v1*3],vert_coords[v1*3+1],vert_coords[v1*3+2]);
		CVector3 vec2 = CVector3(vert_coords[v2*3],vert_coords[v2*3+1],vert_coords[v2*3+2]);
		CVector3 vec3 = CVector3(vert_coords[v3*3],vert_coords[v3*3+1],vert_coords[v3*3+2]);

		CVector3 vec21 = vec2-vec1;
		CVector3 vec31 = vec3-vec1;
		CVector3 cross = vec31.cross(vec21);
		cross.normalize();

		nmap[vec1]+=cross;
		nmap[vec2]+=cross;
		nmap[vec3]+=cross;
	}

	for(int i=0;i<this->vert_count;i++)
	{
		CVector3 v = CVector3(vert_coords[i*3],vert_coords[i*3+1],vert_coords[i*3+2]);
		nmap[v].normalize();
		CVector3 n = nmap[v];
		n=-n;
		vert_norm[i*3] = n.x;
		vert_norm[i*3+1] = n.y;
		vert_norm[i*3+2] = n.z;
	}
	
	//UpdateBTN();
	
}

Thx, but I think that is essentially the same thing that I do, or maybe I didn’t get it. Is there a way to let opengl draw your normals on the surface?

Is there a way to let opengl draw your normals on the surface?

glEnable( GL_AUTO_NORMAL )

Though I have never see it working… You better precompute them manually, it is not that hard.

BTW, your link is broken.

Thx for the fast responses I have implemented the normal drawing thing myself, it wasn’t hard indeed. It looks amazingly cool! I have noticed that where both sides of the terrain come together, they do not match correctly. I’ll try to fix this and then see if this helps.
If you’d have any more suggestions, you’re welcome!

And I have fixed the link.

Ok, I’ve corrected the problems I’ve found, but the main problem is still there. I have uploaded a new picture of it so you can see: http://users.telenet.be/lalaland/wouter/terrein2.jpg
As you can see on this picture, where both ends join, the shading does not match.

Actually I misread your post, you wanted to draw normal vectors and not generate them automatically.

Anyway, I am not sure to understand this:

At the point where the edges of the original surface are brought together you clearly see the edge.

so you mean that you have two surfaces whose one border match exactly and you brought these borders next to each other?
<s>It is quite hard to see on the screenshot, maybe another one taken from another angle may be better.</s>

If I am right, the normal calculation is correct as it appears to be on the rest of the surface. Since normals are interpolated, they are not the same at the surface borders and you get a discontinuity.

Let me explain I start from a .bmp file that contains the information for the height of a 2D surface. Next I want to roll this 2D surface (like creating a cilinder). Thatfore I use cilindric coordinates and join two edges of the 2D surface. Where these edges join is where there is a difference in shading. That’s what you see on the pictures I’ve put online.

Now to have a smooth “cilindric” surface I’ve created the .bmp file so that both edges join smoothly (not only the height values on the edges are identically, also the first derivatives).

So in fact, the difference in shading should not be there…

If you like you can download the code completely from http://users.telenet.be/lalaland/wouter/terrain.rar
Key ‘a’ toggles smooth shading
Key ‘z’ toggles drawing the normals
Key ‘q’ and ‘s’ change the rotate speed
Key ‘w’ and ‘x’ are used to zoom in and out

Ok this is what I suspected.
I did not look your normal computation source yet but I assume that you are computing triangles normals, then intepolating them to calculate vertices normals. Am I correct?

When you interpolate normals at a vertex you must find all triangles that use this vertex then compute the average normal among all triangles that gather this vertex.

The problem is that on the surface edge, vertices are duplicated, therefore, you do not actually find all triangles that “should” gather one vertex.

To solve this problem, the only idea I have is to merge duplicated vertices at the border and change triangle indices to make them use the same vertex at each border side.

GOD!! I made the most stupid error of all mankind!
I thought that I was joining the upper and lower edges of the bmp, but nothing is less true… Rotating my bmp files by 90 degrees solves my problem :s

Thanks anyway for your help and the fast responses!!!
Merging the vertices at the borders is indeed also a good idea to finetune the whole thing. Do you have any idea how to fold the surface once more to map it on a sphere?

Lol! You are welcome. It is strange that you do not see the issue a talked about. Perhaps the surface curvature is too low to notice it.

for the sphere, I have no idea… I would rather generate a sphere mesh and map the heightmap on it to compute the displacement.