Main Page Namespace List Class Hierarchy Compound List File List Namespace Members Compound Members File Members Related Pages

# viewpoints.cpp

    00001 /*==================================================================
    00002
    00003   Program:   View Generation Toolkit
    00004   Module:    viewpoints.cxx
    00005   Language:  C++
    00006   Date:      2002/05/27
    00007   Version:   1.0
    00008
    00009 ===================================================================*/
    00010
    00011 #include <stdio.h>
    00012 #include <stdlib.h>
    00013 #include <math.h>
    00014 #include "viewpoints.h"
    00015
    00016 const double Pi = 3.1415927;
    00017 const double radtodeg = 180.0/Pi;
    00018 const double degtorad = Pi/180.0;
    00019
    00020
    00021 typedef struct {
    00022     double  x, y, z;
    00023 } point;
    00024
    00025 typedef struct {
    00026     point     pt[3];    /* Vertices of triangle */
    00027     double    area;     /* Unused; might be used for adaptive subdivision */
    00028 } triangle;
    00029
    00030 typedef struct {
    00031     int       npoly;    /* # of triangles in object */
    00032     triangle *poly;     /* Triangles */
    00033 } object;
    00034
    00035 typedef struct {
    00036   int    index;
    00037   double value;
    00038 } sortstruct;
    00039
    00040 /* Six equidistant points lying on the unit sphere */
    00041 #define XPLUS {  1,  0,  0 }    /*  X */
    00042 #define XMIN  { -1,  0,  0 }    /* -X */
    00043 #define YPLUS {  0,  1,  0 }    /*  Y */
    00044 #define YMIN  {  0, -1,  0 }    /* -Y */
    00045 #define ZPLUS {  0,  0,  1 }    /*  Z */
    00046 #define ZMIN  {  0,  0, -1 }    /* -Z */
    00047
    00048 /* Vertices of a unit octahedron */
    00049 triangle octahedron[] = {
    00050   { { YPLUS, XPLUS, ZPLUS }, 0.0 },//1
    00051   { { YPLUS, ZPLUS, XMIN  }, 0.0 },//2
    00052   { { YPLUS, XMIN , ZMIN  }, 0.0 },//3
    00053   { { YPLUS, ZMIN , XPLUS }, 0.0 },//4
    00054   { { YMIN , ZPLUS, XPLUS }, 0.0 },//5
    00055   { { YMIN , XMIN , ZPLUS }, 0.0 },//6
    00056   { { YMIN , ZMIN , XMIN  }, 0.0 },//7
    00057   { { YMIN , XPLUS, ZMIN  }, 0.0 } //8
    00058 };
    00059
    00060 /* A unit octahedron */
    00061 object oct = {
    00062     sizeof(octahedron) / sizeof(octahedron[0]),
    00063     &octahedron[0]
    00064 };
    00065
    00066 int compare_sortstruct(const void* c, const void* d)
    00067 {
    00068   sortstruct* a = (sortstruct*)c;
    00069   sortstruct* b = (sortstruct*)d;
    00070
    00071   if(a->value < b->value)
    00072     return -1;
    00073   else if(a->value == b->value)
    00074     return 0;
    00075   else
    00076     return 1;
    00077 }
    00078
    00079 /* Return the midpoint on the line between two points */
    00080 point *midpoint(point* a, point* b)
    00081 {
    00082     static point r;
    00083
    00084     r.x = (a->x + b->x) * 0.5;
    00085     r.y = (a->y + b->y) * 0.5;
    00086     r.z = (a->z + b->z) * 0.5;
    00087
    00088     return &r;
    00089 }
    00090
    00091 double dist(double phi1, double theta1, double phi2, double theta2)
    00092 {
    00093   double cos_alpha;
    00094   phi1 = phi1*degtorad;
    00095   phi2 = phi2*degtorad;
    00096   theta1 = theta1*degtorad;
    00097   theta2 = theta2*degtorad;
    00098
    00099   cos_alpha = sin(phi1)*sin(phi2)*cos(theta1-theta2)+cos(phi1)*cos(phi2);
    00100   return acos(cos_alpha)*radtodeg;
    00101 }
    00102
    00103 point getmidpoint(triangle* t)
    00104 {
    00105         point p;
    00106
    00107         p.x = (t->pt[0].x + t->pt[1].x)/2;
    00108         p.y = (t->pt[0].y + t->pt[1].y)/2;
    00109         p.z = (t->pt[0].z + t->pt[1].z)/2;
    00110
    00111         p.x = (t->pt[2].x - p.x)/3 + p.x;
    00112         p.y = (t->pt[2].y - p.y)/3 + p.y;
    00113         p.z = (t->pt[2].z - p.z)/3 + p.z;
    00114
    00115         return p;
    00116 }
    00117
    00118 /* Normalize a point p */
    00119 point *normalize(point* p)
    00120 {
    00121     static point r;
    00122     double mag;
    00123
    00124     r = *p;
    00125     mag = r.x * r.x + r.y * r.y + r.z * r.z;
    00126     if (mag != 0.0) {
    00127         mag = 1.0 / sqrt(mag);
    00128         r.x *= mag;
    00129         r.y *= mag;
    00130         r.z *= mag;
    00131     }
    00132
    00133     return &r;
    00134 }
    00135
    00136 // Computes r to the power p
    00137 int pow(int r, int p)
    00138 {
    00139   int i = r;
    00140   for(;p > 1;p--) i *= r;
    00141
    00142   return i;
    00143 }
    00144
    00145
    00146 CViewPoints::CViewPoints(int n)
    00147 {
    00148   num = n;
    00149   AllocatePoints(num);
    00150 }
    00151
    00152 CViewPoints::CViewPoints()
    00153 {
    00154     phi = NULL; theta = NULL;
    00155 }
    00156
    00157 CViewPoints::~CViewPoints() {
    00158   // Cleans up memory allocated in ComputePoints
    00159   //
    00160   if(phi != NULL) free(phi);
    00161   if(theta != NULL) free(theta);
    00162 }
    00163
    00164 int* CViewPoints::GetNeighbours(int ind, double range, int& n)
    00165 {
    00166         int i;
    00167         int* index;
    00168         int* retarray = NULL;
    00169
    00170         n = 0;
    00171         if( ( index = (int*)malloc((num-1)*sizeof(int)) ) == NULL) return NULL;
    00172
    00173         for(i=0; i<num; i++) {
    00174                 if(i!=ind && dist(phi[i], theta[i], phi[ind], theta[ind]) <= range) {
    00175                         index[n] = i;
    00176                         n++;
    00177                 }
    00178         }
    00179
    00180         if(n>0 && ( retarray = (int*)malloc(n*sizeof(int)) )!=NULL)
    00181                 for(i=0; i<n; i++) retarray[i] = index[i];
    00182
    00183         free(index);
    00184
    00185         return retarray;
    00186 }
    00187
    00188 int* CViewPoints::GetNeighbours(int ind, int n)
    00189 {
    00190    int i, q;
    00191    sortstruct* ss = NULL;
    00192    int* retarray = NULL;
    00193
    00194    if( (ss = (sortstruct*)malloc((num-1)*sizeof(sortstruct)) ) == NULL)
    00195            return NULL;
    00196
    00197    for(i=0, q=0; i<num; i++) {
    00198      if(i!=ind){
    00199        ss[q].index = i;
    00200        ss[q].value = dist(phi[i], theta[i], phi[ind], theta[ind]);
    00201
    00202        q++;
    00203      }
    00204    }
    00205
    00206    qsort((void *)ss, (size_t)(num-1), sizeof(sortstruct), compare_sortstruct);
    00207
    00208    if( (num-1) < n) n = num-1;
    00209    if( ( retarray = (int*)malloc(n*sizeof(int)) ) == NULL) return NULL;
    00210
    00211    for(i=0; i < n; i++){
    00212      retarray[i] = ss[i].index;
    00213    }
    00214
    00215    free(ss);
    00216
    00217    return retarray;
    00218 }
    00219
    00220 void CViewPoints::AllocatePoints(int n)
    00221 {
    00222   if( (phi = (float*)malloc(sizeof(float)*n)) == NULL) Abort();
    00223   if( (theta = (float*)malloc(sizeof(float)*n)) == NULL){
    00224     free(phi);
    00225     Abort();
    00226   }
    00227 }
    00228
    00229 void CViewPoints::Abort()
    00230 {
    00231   // Cleans up memory allocated in ComputePoints
    00232   //
    00233   if(phi != NULL) free(phi);
    00234   if(theta != NULL) free(theta);
    00235
    00236   printf("Abort has occured\n");
    00237
    00238   abort();
    00239 }
    00240
    00241 float CViewPoints::Phi(int ind) {
    00242   if((ind > -1) && (ind < num))
    00243         return phi[ind];
    00244   else
    00245         return -1;
    00246 }
    00247
    00248 float CViewPoints::Theta(int ind) {
    00249   if((ind > -1) && (ind < num))
    00250           return theta[ind];
    00251   else
    00252           return -1;
    00253 }
    00254
    00255 int CViewPoints::Size()
    00256 {
    00257   return num;
    00258 }
    00259
    00260 int CViewPoints::Level()
    00261 {
    00262   return mLevel;
    00263 }
    00264
    00265 COctahedronPoints::COctahedronPoints(int n)
    00266 {
    00267   mLevel = n;
    00268
    00269   if(n < 1) abort();
    00270
    00271   num = 2*pow(4,n);
    00272
    00273   AllocatePoints(num);
    00274   ComputePoints();
    00275 }
    00276
    00277
    00278 void COctahedronPoints::ComputePoints()
    00279 {
    00280     object *old = &oct,         /* Default is octahedron */
    00281            *newm;
    00282
    00283     int     i,
    00284             level,              /* Current subdivision level */
    00285             maxlevel = mLevel;  /* Maximum subdivision level */
    00286
    00287
    00288     /* Subdivide each starting triangle (maxlevel - 1) times */
    00289     for (level = 1; level < maxlevel; level++) {
    00290         /* Allocate a new object */
    00291         newm = (object *)malloc(sizeof(object));
    00292         if (newm == NULL) {
    00293             Abort();
    00294         }
    00295         newm->npoly = old->npoly * 4;
    00296
    00297         /* Allocate 4* the number of points in the current approximation */
    00298         newm->poly  = (triangle *)malloc(newm->npoly * sizeof(triangle));
    00299         if (newm->poly == NULL) {
    00300             Abort();
    00301         }
    00302
    00303         /* Subdivide each triangle in the old approximation and normalize
    00304          *  the new points thus generated to lie on the surface of the unit
    00305          *  sphere.
    00306          * Each input triangle with vertices labelled [0,1,2] as shown
    00307          *  below will be turned into four new triangles:
    00308          *
    00309          *                      Make new points
    00310          *                          a = (0+2)/2
    00311          *                          b = (0+1)/2
    00312          *                          c = (1+2)/2
    00313          *        1
    00314          *       /\             Normalize a, b, c
    00315          *      /  \
    00316          *    b/____\ c         Construct new triangles
    00317          *    /\    /\              [0,b,a]
    00318          *   /  \  /  \             [b,1,c]
    00319          *  /____\/____\            [a,b,c]
    00320          * 0      a     2           [a,c,2]
    00321          */
    00322         for (i = 0; i < old->npoly; i++) {
    00323             triangle
    00324                  *oldt = &old->poly[i],
    00325                  *newt = &newm->poly[i*4];
    00326             point a, b, c;
    00327
    00328             a = *normalize(midpoint(&oldt->pt[0], &oldt->pt[1]));
    00329             b = *normalize(midpoint(&oldt->pt[0], &oldt->pt[2]));
    00330             c = *normalize(midpoint(&oldt->pt[1], &oldt->pt[2]));
    00331
    00332             newt->pt[0] = oldt->pt[0];
    00333             newt->pt[1] = a;
    00334             newt->pt[2] = b;
    00335             newt++;
    00336
    00337             newt->pt[0] = a;
    00338             newt->pt[1] = oldt->pt[1];
    00339             newt->pt[2] = c;
    00340             newt++;
    00341
    00342             newt->pt[0] = c;
    00343             newt->pt[1] = a;
    00344             newt->pt[2] = b;
    00345             newt++;
    00346
    00347             newt->pt[0] = b;
    00348             newt->pt[1] = c;
    00349             newt->pt[2] = oldt->pt[2];
    00350         }
    00351
    00352         if (level > 1) {
    00353             free(old->poly);
    00354             free(old);
    00355         }
    00356
    00357         /* Continue subdividing new triangles */
    00358         old = newm;
    00359     }
    00360
    00361
    00362     /* Spit out coordinates for each triangle */
    00363     point p;
    00364     double l1, l2;
    00365
    00366     triangle *t;
    00367     for (i = 0; i < old->npoly; i++) {
    00368         t = &old->poly[i];
    00369
    00370         p = getmidpoint(t);
    00371         l1 = sqrt(p.x*p.x + p.y*p.y + p.z*p.z);
    00372         l2 = sqrt(p.x*p.x + p.y*p.y);
    00373
    00374         phi[i] = acos(p.z/l1)*radtodeg;
    00375         if(l2 == 0)
    00376           theta[i] = 0;
    00377         else {
    00378           theta[i] = acos(p.x/l2)*radtodeg;
    00379           if(p.y < 0) theta[i] = 360 - theta[i];
    00380         }
    00381     }
    00382
    00383     free(old);
    00384 }