vtkRecursiveSphereDirectionEncoder.cxx 11.9 KB
Newer Older
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkRecursiveSphereDirectionEncoder.cxx

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9

10 11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
13 14 15

=========================================================================*/
#include "vtkRecursiveSphereDirectionEncoder.h"
16 17
#include "vtkObjectFactory.h"

18
#include <cmath>
19

Brad King's avatar
Brad King committed
20
vtkStandardNewMacro(vtkRecursiveSphereDirectionEncoder);
21 22 23 24 25 26 27

// Construct the object. Initialize the index table which will be
// used to map the normal into a patch on the recursively subdivided
// sphere.
vtkRecursiveSphereDirectionEncoder::vtkRecursiveSphereDirectionEncoder()
{
  this->RecursionDepth = 6;
28 29
  this->IndexTable = nullptr;
  this->DecodedNormal = nullptr;
30 31 32 33 34 35
  this->InitializeIndexTable();
}

// Destruct a vtkRecursiveSphereDirectionEncoder - free up any memory used
vtkRecursiveSphereDirectionEncoder::~vtkRecursiveSphereDirectionEncoder()
{
36 37
  delete [] this->IndexTable;
  delete [] this->DecodedNormal;
38 39 40 41 42 43 44 45
}


int vtkRecursiveSphereDirectionEncoder::GetEncodedDirection( float n[3] )
{
  float t;
  int   value;
  int   xindex, yindex;
46
  float x, y;
47 48

  if ( this->IndexTableRecursionDepth != this->RecursionDepth )
49
  {
50
    this->InitializeIndexTable();
51
  }
52 53

  // Convert the gradient direction into an encoded index value
54
  // This is done by computing the (x,y) grid position of this
55 56
  // normal in the 2*NORM_SQR_SIZE - 1 grid, then passing this
  // through the IndexTable to look up the 16 bit index value
57

58 59
  // Don't use fabs because it is slow - just convert to absolute
  // using a simple conditional.
60 61 62
  t =
    ((n[0]>=0.0)?(n[0]):(-n[0])) +
    ((n[1]>=0.0)?(n[1]):(-n[1])) +
63
    ((n[2]>=0.0)?(n[2]):(-n[2]));
64

65
  if ( t )
66
  {
67

68
    t = 1.0 / t;
69

70 71
    x = n[0] * t;
    y = n[1] * t;
72

73
    xindex = (int)((x+1.0)*(float)(this->InnerSize) + 0.5);
74
    yindex = (int)((y+1.0)*(float)(this->InnerSize) + 0.5);
75

76
    if ( xindex > 2*this->InnerSize )
77
    {
78
      xindex = 2*this->InnerSize;
79
    }
80
    if ( yindex > 2*this->InnerSize )
81
    {
82
      yindex = 2*this->InnerSize;
83
    }
84

85
    value = this->IndexTable[xindex*(this->OuterSize+this->InnerSize) + yindex];
86

87
    // If the z component is less than 0.0, add this->GridSize to the
88
    // index
Bill Lorensen's avatar
Bill Lorensen committed
89
    if ( n[2] < 0.0 )
90
    {
91
      value += this->GridSize;
92
    }
93
  }
94
  else
95
  {
96
    value = 2*this->GridSize;
97
  }
98 99 100

  return value;
}
101

102 103 104
float *vtkRecursiveSphereDirectionEncoder::GetDecodedGradient( int value )
{
  if ( this->IndexTableRecursionDepth != this->RecursionDepth )
105
  {
106
    this->InitializeIndexTable();
107
  }
108 109 110 111

  return (this->DecodedNormal + value*3);
}

112
int vtkRecursiveSphereDirectionEncoder::GetNumberOfEncodedDirections( )
113 114 115 116
{
  int     outer_size, inner_size;
  int     norm_size;

Jim Miller's avatar
Jim Miller committed
117
  outer_size = (int)(pow( 2.0, (double) this->RecursionDepth ) + 1);
118 119 120 121 122 123 124
  inner_size = outer_size - 1;

  norm_size = outer_size * outer_size + inner_size * inner_size;

  return (norm_size*2 + 1);
}

125
float *vtkRecursiveSphereDirectionEncoder::GetDecodedGradientTable( )
126 127
{
  if ( this->IndexTableRecursionDepth != this->RecursionDepth )
128
  {
129
    this->InitializeIndexTable();
130
  }
131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147

  return this->DecodedNormal;
}

// Initialize the index table.  This is a 2*NORM_SQR_SIZE - 1 by
// 2*NORM_SQR_SIZE - 1 entry table that maps (x,y) grid position to
// encoded normal index.  The grid position is obtained by starting
// with an octahedron (comprised of 8 triangles forming a double
// pyramid). Each triangle is then replaced by 4 triangles by joining
// edge midpoints.  This is done recursively until NORM_SQR_SIZE
// vertices exist on each original edge. If you "squish" this octahedron,
// it will look like a diamond.  Then rotate it 45 degrees, it will
// look like a square.  Then look at the pattern of vertices - there
// is a NORM_SQR_SIZE by NORM_SQR_SIZE grid, with a (NORM_SQR_SIZE-1) by
// NORM_SQR_SIZE - 1 grid inside of it.  The vertices all fall on
// (x,y) locatiions in a grid that is 2*NORM_SQR_SIZE - 1 by
// 2*NORM_SQR_SIZE - 1, although not every (x,y) location has a vertex.
148
void vtkRecursiveSphereDirectionEncoder::InitializeIndexTable( )
149 150 151 152 153
{
  int     i, j, index, max_index;
  int     xindex, yindex;
  float   x, y, z, tmp_x, tmp_y;
  float   norm;
154
  int     limit;
155 156

  // Free up any memory previously used
157 158
  delete [] this->IndexTable;
  delete [] this->DecodedNormal;
159

Jim Miller's avatar
Jim Miller committed
160
  this->OuterSize = (int)(pow( 2.0, (double) this->RecursionDepth ) + 1);
161
  this->InnerSize = this->OuterSize - 1;
162
  this->GridSize =
163 164 165
    this->OuterSize * this->OuterSize +
    this->InnerSize * this->InnerSize;

166 167

  // Create space for the tables
168
  this->IndexTable = new int [(this->OuterSize + this->InnerSize) *
169
                              (this->OuterSize + this->InnerSize)];
170 171 172

  // Initialize the table to -1 -- we'll use this later to determine which
  // entries are still not filled in
173
  for ( i = 0; i < ( (this->OuterSize + this->InnerSize) *
174
                     (this->OuterSize + this->InnerSize) ); i ++ )
175
  {
176
      this->IndexTable[i] = -1;
177
  }
178

179 180
  this->DecodedNormal =
    new float [ 3 * ( 1 +
181 182
                      2 * this->OuterSize * this->OuterSize +
                      2 * this->InnerSize * this->InnerSize ) ];
183 184 185 186 187 188 189 190 191 192

  // Initialize the index
  index = 0;

  // max_index indicates the largest index we will get - the number
  // of vertices in the two-grid square.  This represents half the
  // normals, and max_index is used to offset from one half into the
  // other.  One half of the normals have z components >= 0, and the
  // second half (all with indices above max_index) have z components
  // that are <= 0.
193
  max_index =  this->GridSize;
194 195 196 197 198 199

  // The last normal (max_index*2) is the zero normal
  this->DecodedNormal[3*(max_index*2) + 0] = 0.0;
  this->DecodedNormal[3*(max_index*2) + 1] = 0.0;
  this->DecodedNormal[3*(max_index*2) + 2] = 0.0;

200 201
  // The outer loop is for this->OuterSize + this->InnerSize rows
  for ( i = 0; i < this->OuterSize + this->InnerSize; i++ )
202
  {
203
    // Compute the y component for this row
204
    tmp_y = (float)(2*i)/(float)(this->InnerSize*2) - 1.0;
205 206

    // On the odd rows, we are doing the small grid which has
207 208
    // this->InnerSize elements in it
    limit = ( i%2 )?(this->InnerSize):(this->OuterSize);
209 210

    for ( j = 0; j < limit; j++ )
211
    {
212 213
      // compute the x component for this column
      if ( i%2 )
214
      {
215
        tmp_x = (float)(2*j)/(float)(this->InnerSize) -
216
          1.0 + (1.0/(float)(this->InnerSize));
217
      }
218
      else
219
      {
220
        tmp_x = (float)(2*j)/(float)(this->InnerSize) - 1.0;
221
      }
222

223
      // rotate by 45 degrees
224
      // This rotation intentionally does not preserve length -
225 226 227 228
      // we could have tmp_x = 1.0 and tmp_y = 1.0, we want this
      // to lie within [-1.0,1.0] after rotation.
      x = 0.5 * tmp_x - 0.5 * tmp_y;
      y = 0.5 * tmp_x + 0.5 * tmp_y;
229

230 231
      // compute the z based on the x and y values
      if ( x >= 0 && y >= 0 )
232
      {
233
        z = 1.0 - x - y;
234
      }
235
      else if ( x >= 0 && y < 0 )
236
      {
237
        z = 1.0 - x + y;
238
      }
239
      else if ( x < 0 && y < 0 )
240
      {
241
        z = 1.0 + x + y;
242
      }
243
      else
244
      {
245
        z = 1.0 + x - y;
246
      }
247

248
      // Normalize this direction and set the DecodedNormal table for
249
      // this index to this normal.  Also set the corresponding
250 251 252 253 254 255 256 257 258 259 260
      // entry for this normal with a negative z component
      norm = sqrt( (double)( x*x + y*y + z*z ) );
      this->DecodedNormal[3*index + 0] = x / norm;
      this->DecodedNormal[3*index + 1] = y / norm;
      this->DecodedNormal[3*index + 2] = z / norm;
      this->DecodedNormal[3*(index+max_index) + 0] =   x / norm;
      this->DecodedNormal[3*(index+max_index) + 1] =   y / norm;
      this->DecodedNormal[3*(index+max_index) + 2] = -(z / norm);

      // Figure out the location in the IndexTable. Be careful with
      // boundary conditions.
261
      xindex = (int)((x+1.0)*(float)(this->InnerSize) + 0.5);
262 263
      yindex = (int)((y+1.0)*(float)(this->InnerSize) + 0.5);
      if ( xindex > 2*this->InnerSize )
264
      {
265
        xindex = 2*this->InnerSize;
266
      }
267
      if ( yindex > 2*this->InnerSize )
268
      {
269
        yindex = 2*this->InnerSize;
270
      }
271
      this->IndexTable[xindex*(this->OuterSize+this->InnerSize) + yindex] = index;
272 273 274 275 276 277 278 279 280 281

      // Do the grid location to the left - unless we are at the left
      // border of the grid. We are computing indices only for the
      // actual vertices of the subdivided octahedron, but we'll
      // convert these into the IndexTable coordinates and fill in
      // the index for the intermediate points on the grid as well.
      // This way we can't get bitten by a scan-conversion problem
      // where we skip over some table index due to precision, and
      // therefore it doesn't have a valid value in it.
      if ( tmp_x > 0 )
282
      {
283 284
        x = 0.5 * (tmp_x - (1.0/(float)this->InnerSize)) - 0.5 * tmp_y;
        y = 0.5 * (tmp_x - (1.0/(float)this->InnerSize)) + 0.5 * tmp_y;
285
        xindex = (int)((x+1.0)*(float)(this->InnerSize) + 0.5);
286 287
        yindex = (int)((y+1.0)*(float)(this->InnerSize) + 0.5);
        if ( xindex > 2*this->InnerSize )
288
        {
289
          xindex = 2*this->InnerSize;
290
        }
291
        if ( yindex > 2*this->InnerSize )
292
        {
293 294
          yindex = 2*this->InnerSize;
        }
295 296
        this->IndexTable[xindex*(this->OuterSize+this->InnerSize) + yindex] = index;
      }
297

298 299 300
      // On the odd rows we also need to do the last grid location on
      // the right.
      if ( (i%2) && (j == limit - 1) )
301
      {
302 303
        x = 0.5 * (tmp_x + (1.0/(float)this->InnerSize)) - 0.5 * tmp_y;
        y = 0.5 * (tmp_x + (1.0/(float)this->InnerSize)) + 0.5 * tmp_y;
304
        xindex = (int)((x+1.0)*(float)(this->InnerSize) + 0.5);
305 306
        yindex = (int)((y+1.0)*(float)(this->InnerSize) + 0.5);
        if ( xindex > 2*this->InnerSize )
307
        {
308
          xindex = 2*this->InnerSize;
309
        }
310
        if ( yindex > 2*this->InnerSize )
311
        {
312 313
          yindex = 2*this->InnerSize;
        }
314 315
        this->IndexTable[xindex*(this->OuterSize+this->InnerSize) + yindex] = index;
      }
316

317 318
      // Increment the index
      index++;
319
    }
320
  }
321

322

323 324 325
  // The index table has been initialized for the current recursion
  // depth
  this->IndexTableRecursionDepth = this->RecursionDepth;
326 327 328 329 330 331 332


  // Spread the first index value in each row to the left, and the last to the right.
  // This is because we have only filled in a diamond of index values within the square
  // grid, and we need to be careful at the edges due to precision problems. This way
  // we won't be able to access a table location that does not have a valid index in it.
  for ( j = 0; j < this->OuterSize + this->InnerSize; j++ )
333
  {
334 335
    // Start from the middle going right, copy the value from the left if
    // this entry is not initialized
336
    for ( i = (this->OuterSize+this->InnerSize)/2;
337
          i < this->OuterSize + this->InnerSize; i++ )
338
    {
339
      if ( this->IndexTable[j*(this->OuterSize+this->InnerSize)+i] == -1 )
340
      {
341 342
        this->IndexTable[j*(this->OuterSize+this->InnerSize)+i] =
          this->IndexTable[j*(this->OuterSize+this->InnerSize)+i-1];
343
      }
344
    }
345 346 347 348

    // Start from the middle going left, copy the value from the right if
    // this entry is not initialized
    for ( i = (this->OuterSize+this->InnerSize)/2; i >= 0; i-- )
349
    {
350
      if ( this->IndexTable[j*(this->OuterSize+this->InnerSize)+i] == -1 )
351
      {
352 353
        this->IndexTable[j*(this->OuterSize+this->InnerSize)+i] =
          this->IndexTable[j*(this->OuterSize+this->InnerSize)+i+1];
354 355
      }
    }
356
  }
357 358 359 360
}


// Print the vtkRecursiveSphereDirectionEncoder
361
void vtkRecursiveSphereDirectionEncoder::PrintSelf(ostream& os, vtkIndent indent)
362
{
Brad King's avatar
Brad King committed
363
  this->Superclass::PrintSelf(os,indent);
364

365
  os << indent << "Number of encoded directions: " <<
366
    this->GetNumberOfEncodedDirections() << endl;
367

368
  os << indent << "Recursion depth: " << this->RecursionDepth << endl;
369
}