# Easter 2023: Superquadric and Squared-Eggs

First there was the Great Cosmic Egg.

Huai-nan Tzu, China 100 BCE

Another Easter is arriving, and here I am for a new article on old friends: the Eggs. In my blog for Easter 2021, I mentioned squared eggs paraphrasing a sentence attributed to the extraordinary jeweler Peter Carl Fabergé:” This year, your Highness, we will be featuring square eggs.” A hen is unlikely to make a cubic-shaped egg, but we can still transform a cooked egg to get the cuboid shape. If you look for squared eggs on the internet, you will find a tool called the egg cuber that does the magic. Unfortunately, I could not yet try one, but it seems to be doing an excellent job from the reviews. The shape of the egg is not a perfect cube but a cuboid, namely a cube with rounded corners. In my previous article, I anticipated that this shape could be obtained using the superquadric function. Curiouser and curiouser! In this article, I will give a bit of mathematical background and even a source code to play with such a function as my Easter Bunny gift.

So what is superquadric function? A superellipsoid is a generic definition for a large class of solids with Lamé curves (superellipses) as section profile. Gabriel Lamé (1795-1870) was a French mathematician and engineer who studied superellipses, two-dimensional curves resembling deformed ellipses and hyperbolas . In the Figure 1., the implicit functions and some example of these curves are reported

By rotating these curves, generating the so-called Lamé surfaces, the superellipsoid mentioned at the beginning is possible. This function is a graphical primitive, controlled by two parameters representing a sphere, square box, or closed cylindrical containers [3,4].

The implicit equations of a super ellipsoid is $F(x,y,z)=\left [\left( \frac{x}{a_x}\right)^{2/\eta_2}+\left( \frac{y}{a_y}\right)^{2/\eta_2} \right]^{\eta_2/\eta_1} + \left( \frac{z}{a_z}\right)^{2/\eta_2}$

with $a_x,a_y,a_z$, the semi-diameters of the resulting solid and $\eta_1$ and $\eta_2$  are two positive real numbers used as parameters to control the squareness of the shape.

The parametric representation of the curve is given by the following equations $x(\theta,\omega) =a_x \cos(\theta)^{\eta_1}\cos(\omega)^{\eta_2}$ $y(\theta,\omega) =a_y \cos(\theta)^{\eta_1}\sin(\omega)^{\eta_2}$ $z(\theta,\omega) =a_z \sin(\theta)^{\eta_1}$

with the angle parameters are defined as $\theta \in [-\pi/2,\pi/2]$ and $\omega \in [-\pi,\pi].$

The last equation can be easily implemented in a C++ graphical program using the OpenGL graphical library. The source code in the appendix is my Easter gift! I have reused the codes presented in my previous articles on the egg shape representation.

And now finally we can show how a rainbowly square egg look like using $\eta_1=0.5$ and $\eta_2=0.5$ !

## REFERENCES

1. G. Lamé,. Examen des différentes méthodes employées pour résoudre les problèmes de géométrie, 1918, Vve Courcier (paris). (see at https://gallica.bnf.fr/ark:/12148/bpt6k92728m/f139.item#)
2. See https://mathcurve.com/courbes2d.gb/lame/lame.shtml
3. Alan H. Barr. Superquadrics and Angle-Preserving Transformations. (1981) IEEE Computer Graphics and Applications1(1): 11–23.
4. Alan H. Barr. Rigid Physically Based Superquadrics. pp. 137–159 in  Graphics Gems III, edited by D. Kirk.
5. D. Hearn and M.P. Backer. Computer Graphics with OpenGL. Third Edition. 2004, Pearson.

## APPENDIX

Superellipsoid curve generator. The program is written in C++ and compiled using gnu g++ compiler on a MacBook Pro 2021 (Apple M1 Pro processor) under macOS Ventura. Please read the licence and disclaimer in the About page of this blog about the program provided in this site.

Command used for the compilation.

g++ -arch arm64 SquareEggs.cpp -framework opengl -framework glut -Wno-deprecated -L../glui-2.36/src/lib/ -lglui -o SquareEggs.exe

Note that I have installed the glui library (you can download it from github) in a local folder.

SOURCE CODE

#include <fstream>
#include <iostream>
#include <iomanip>
#include <math.h>
#include <GLUT/glut.h>
#include <OpenGL/gl.h>
#include <OpenGL/glu.h>
#include <GL/glui.h>

#define PI 3.14159265358979323846
#define TWOPI 2*3.14159265358979323846
#define PID2 3.14159265358979323846/2.0
#define sign(a) ( ( (a) < 0 )  ?  -1   : ( (a) > 0 ) )

// Define superellipsoid parameters

float rx = 1.0;   // x-radius
float ry = 1.0;   // y-radius
float rz = 1.0;   // z-radius
float n1 = 2.0;  // first exponent
float n2 = 2.0;  // second exponent

// Define resolution
int res = 30;

int   main_window;
GLint spin= 0;
GLint spin1=0;
GLint egg= 0;
GLUI_Spinner    *spi;

//Viewer options (GluLookAt)

float fovy = 60.0, aspect = 1.0, zNear = 1.0, zFar = 100.0;

//Mouse modifiers

float depth = -1;
float scale=2.5;
float psi=0, theta=0;
float downX, downY;
bool leftButton = false, middleButton = false, rightButton=false;
int  showAxis;
int  showLight=1,showLight1=1;
int  sqstripes=-1;
int  chir;

void display(void);
int Squaredegg(void);

void myGlutIdle( void )
{
if ( glutGetWindow() != main_window )
glutSetWindow(main_window);

glutPostRedisplay();
}
/*
============ 3D math routines
*/

//------ Returns the cross product of 2 vectors

void CrossProduct (double M[], double N[], double CS[]) {
CS=M * N - M * N;
CS=M * N - M * N;
CS=M * N - M * N;
}

//------ Returns the length of a vector

double GetVectorLength (double M[]) {
return sqrt( M * M + M * M + M * M );
}

//------ Returns the sum of two vectors

void   AddVectors (double M[], double N[], double R[]){
R= M + N;
R= M + N;
R= M + N;
}

//------ Returns the vector scaled by the last parameter

void ScaleVector (double M[], double a, double N[]) {
N= M * a;
N= M * a;
N= M * a;
}

//------ Returns a normalized vector (length = 1)

void NormalizeVector (double M[],double R) {
double norm = GetVectorLength(M);
if (norm == 0) norm  =1.0;
ScaleVector( M, 1./ norm, R );
}

//------ Returns the unit normal vector of a triangle specified by the three
// points P1, P2, and P3.

void FindUnitNormal (double P1,double P2,double P3,double N){
double D1,D2;
double R;
D1=P1-P2;
D1=P1-P2;
D1=P1-P2;
D2=P2-P3;
D2=P2-P3;
D2=P2-P3;
CrossProduct(D1,D2,R);
NormalizeVector(R,N);
}

//------ Initialization routine

void init () {

GLfloat  mat_ambient[]  = { 0.25,     0.20725,  0.20725,  0 };
GLfloat  mat_diffuse[]  = { 1.0,        0.829,    0.829,    0 };
GLfloat  mat_specular[] = { 0.296648, 0.296648, 0.296648, 0 };
GLfloat  light_diffuse[]  = { 0.8, 0.8, 0.8, 1.0 };
GLfloat  light_ambient[]  = { 0.8, 0.8, 0.8, 1.0 };
GLfloat  light_specular[] = { 0.5, 0.5, 0.5, 1.0 };

GLfloat  light0_position[] = { -20.0, 20.0, 0.0, 0 };
GLfloat  light1_position[] = { 10.0, 50.0, 0.0, 0 };

GLfloat  shininess=0.088 * 128;

glClearColor( 0  , 0  , 0  , 0   );    // Black background
glEnable(GL_MULTISAMPLE);      // Enable multisample antialiasing
glEnable(GL_DEPTH_TEST);       // Enable hidden surface removal
glEnable(GL_LIGHTING);

// set light 0

glEnable(GL_LIGHT0);
glLightfv( GL_LIGHT0, GL_POSITION, light0_position );
glLightfv( GL_LIGHT0, GL_DIFFUSE,  light_diffuse );
glLightfv( GL_LIGHT0, GL_AMBIENT,  light_ambient );
glLightfv( GL_LIGHT0, GL_SPECULAR, light_specular );

// set light 1

glEnable(GL_LIGHT1);
glLightfv( GL_LIGHT1, GL_POSITION, light1_position );
glLightfv( GL_LIGHT1, GL_DIFFUSE,  light_diffuse );
glLightfv( GL_LIGHT1, GL_AMBIENT,  light_ambient );
glLightfv( GL_LIGHT1, GL_SPECULAR, light_specular );

// set material

glMaterialfv( GL_FRONT, GL_AMBIENT,  mat_ambient );
glMaterialfv( GL_FRONT, GL_DIFFUSE,  mat_diffuse );
glMaterialfv( GL_FRONT, GL_SPECULAR, light_specular );
glMaterialf( GL_FRONT, GL_SHININESS, shininess );

egg=Squaredegg();
}

void display () {

// Clear window

glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

gluLookAt( 2, 4, 10, 0, 0, 0, 0, 1, 0 );
glPushMatrix();
glScalef( scale, scale, scale );

//Motion Options

glTranslatef(0.0, 0.0, -depth);
glRotatef(-theta, 1.0, 0.0, 0.0);
glRotatef(psi, 0.0, 1.0, 0.0);

if (showLight) {
glEnable(GL_LIGHT0);
}
else
{
glDisable(GL_LIGHT0);
}

if (showLight1) {
glEnable(GL_LIGHT1);
}
else
{
glDisable(GL_LIGHT1);
}

glCallList(egg);            // draw the squared eggs

glPopMatrix();
glutSwapBuffers();
}

void SuperEllipsoidCurve(double theta1,double theta2,double p1,double p2, double x[])
{
double tmp;
double ct1,ct2,st1,st2;

ct1 = cos(theta1);
ct2 = cos(theta2);
st1 = sin(theta1);
st2 = sin(theta2);

tmp  = sign(ct1) * pow(fabs(ct1),p1);
x = tmp * sign(ct2) * pow(fabs(ct2),p2);
x = sign(st1) * pow(fabs(st1),p1);
x = tmp * sign(st2) * pow(fabs(st2),p2);
}

int Squaredegg(void)
{
GLuint  egg;

double x,x1,x2,NN;
double theta3;

// Draw superellipsoid

egg   = glGenLists(1);
glNewList( egg, GL_COMPILE );

double dt = 0.01 * TWOPI / res;
for (int j=0;j<res/2;j++) {
double theta1 = j * TWOPI / (double)res - PID2;
double theta2 = (j + 1) * TWOPI / (double)res - PID2;

for (int i=0;i<=res;i++) {
double r1= rand() / double(RAND_MAX);
double r2= rand() / double(RAND_MAX);
double r3= rand() / double(RAND_MAX);

GLfloat  mat_ambient[]  = { r1 ,     r2,  r3,  0 };
glMaterialfv( GL_FRONT, GL_AMBIENT,  mat_ambient );

if (i == 0 || i == res)
theta3 = 0;
else
theta3 = i * TWOPI / res;

SuperEllipsoidCurve(theta2,theta3,n1,n2, x);
SuperEllipsoidCurve(theta2+dt,theta3,n1,n2, x1);
SuperEllipsoidCurve(theta2,theta3+dt,n1,n2, x2);
FindUnitNormal (x1,x,x2,NN);
glNormal3d(NN,NN,NN);

glTexCoord2f(i/(double)res,2*(j+1)/(double)res);
glVertex3f(rx*x,ry*x,rz*x);

SuperEllipsoidCurve(theta1,theta3,n1,n2, x);
SuperEllipsoidCurve(theta1+dt,theta3,n1,n2, x1);
SuperEllipsoidCurve(theta1,theta3+dt,n1,n2, x2);
FindUnitNormal (x1,x,x2,NN);
glNormal3d(NN,NN,NN);

glTexCoord2f(i/(double)res,2*j/(double)res);
glVertex3f(rx*x,ry*x,rz*x);

}
glEnd();
}

glEndList();

return egg;

}

void reshape (int w, int h)
{
if (h == 0 || w == 0) return;

glViewport( 0, 0, w, h );
glMatrixMode(GL_PROJECTION);
gluPerspective( 60, h ? w / h : 0, 1, 20 );
glMatrixMode(GL_MODELVIEW);
}

//------ Routine for rotating the scene

void spinDisplay (){

int TimeNow = glutGet(GLUT_ELAPSED_TIME);
int WaitUntil=0;

if ( TimeNow >= WaitUntil ) {
spin += 1;

if ( spin > 360 )spin = spin - 360;
glutPostRedisplay();
WaitUntil = TimeNow + 1000 / 25;    // 25 frames/s
}
}

/* Callbacks */
void mouseCallback(int button, int state, int x, int y)
{
downX = x; downY = y;
leftButton = ((button == GLUT_LEFT_BUTTON) && (state == GLUT_DOWN));
rightButton = ((button == GLUT_RIGHT_BUTTON) && (state == GLUT_DOWN));
middleButton = ((button == GLUT_MIDDLE_BUTTON) &&  (state == GLUT_DOWN));
}

void motionCallback(int x, int y)
{
if (leftButton) //Rotate
{
psi += (x-downX)/4.0;
theta += (downY-y)/4.0;
}
if (rightButton) //Scale
{
if (depth + (downY - y)/10.0 < zFar-10 && depth + (downY - y)/10.0 > zNear+3)
depth += (downY - y)/10.0;
}
downX = x;
downY = y;

glutPostRedisplay();
}

void control_cb( int control )
{
egg=Squaredegg();
}

int main (int argc, char **argv)
{
//Initialize GLUT
glutInit(&argc, argv);
//double buffering used to avoid flickering problem in animation
glutInitDisplayMode(
GLUT_DOUBLE         // Double buffering
| GLUT_RGB            // RGB color mode
| GLUT_DEPTH          // Hidden surface removal
| GLUT_MULTISAMPLE    // Multisample antialiasing
);

// window size
glutInitWindowSize( 600, 600 );
// create the window
main_window = glutCreateWindow("Egg shapes");
init();
//Assign  the function used in events
glutDisplayFunc(display);

glutReshapeFunc( reshape );
// mouse motion
glutMouseFunc(mouseCallback);
glutMotionFunc(motionCallback);

/****************************************/

GLUI *glui = GLUI_Master.create_glui( "Command & Control",300,500 );

// Define resolution
int res = 30;
spi= new GLUI_Spinner( glui, "X-radius (rx):", &rx,200,control_cb);
spi ->set_float_limits( 0.1, 10. );
spi->set_speed(0.1);
spi= new GLUI_Spinner( glui, "Y-radius (ry):", &ry,200,control_cb);
spi ->set_float_limits( 0.1, 10. );
spi->set_speed(0.1);
spi= new GLUI_Spinner( glui, "Z-radius (rz):", &rz,200,control_cb);
spi ->set_float_limits( 0.1, 10. );
spi->set_speed(0.1);
spi= new GLUI_Spinner( glui, "First exponent (n1):", &n1,200,control_cb);
spi ->set_float_limits( 0.1, 10. );
spi->set_speed(0.1);
spi= new GLUI_Spinner( glui, "Second exponent (n2):", &n2,200,control_cb);
spi ->set_float_limits( 0.1, 10. );
spi->set_speed(0.05);
spi= new GLUI_Spinner( glui, "Resolution (res):", &res,200,control_cb);
spi ->set_float_limits( 30, 300. );
spi->set_speed(0.1);

new GLUI_Checkbox( glui, "Show Lights 0     ", &showLight);
new GLUI_Checkbox( glui, "Show Lights 1    ", &showLight1);

(new GLUI_Spinner( glui, "Scaling:", &scale,206, control_cb )) ->set_float_limits( 0, 20. );

glui->set_main_gfx_window( main_window );

/* We register the idle callback with GLUI, *not* with GLUT */
GLUI_Master.set_glutIdleFunc( myGlutIdle );
//Let start glut loop
glutMainLoop();
return 0;
}



This site uses Akismet to reduce spam. Learn how your comment data is processed.