//...................................................................
//
// surface.cpp
// INCLUDES
#include <stdio.h> // For IO
#include <math.h> // For Math
#include "surface.h" // For Surface class
// DEFINES
#define IMAX_INTENSITY 255 // Maximum display intensity
#define RAD 0.01745329 // Degrees to radian conversion
#define ILLUM_LAT 45.0 // Illumination elevation angle
#define ILLUM_LON 45.0 // Illumination azimuth angle
//.............................. CSurface ...........................
// Constructor
//
CSurface::CSurface(void)
{
m_GridSize = m_Smax = 0; // zero class members
m_Nx = m_Ny = m_Nz = m_Nb = 0;
m_Image = 0;
m_Hz = 0;
m_Elevations = 0;
}
//............................. ~CSurface ...........................
// Destructor
//
CSurface::~CSurface(void)
{
m_GridSize = m_Smax = 0;
m_Nx = m_Ny = m_Nz = m_Nb = 0; // zero class members
if(m_Hz != 0) delete [] m_Hz; // free alloc'd memory
if(m_Image != 0) delete [] m_Image;
if(m_Elevations != 0) delete [] m_Elevations;
m_Image = 0;
m_Hz = 0;
m_Elevations = 0;
}
//.......................... surfaceRead ............................
// Read elevation data.
//
short CSurface::surfaceRead(char *fileName)
{
// Application dependent. Include code to read surface data. This
// function should, at a minimum, allocate memory for m_Elevations
// and fill this memory with data defining the surface. The
// surface dimensions m_Nx and m_Ny should also be set here along
// with m_GridSize the uniform spacing between elevation data,
// and the bitmap edge buffer m_Nb. If no edge buffer is required
// m_Nb should be set to one.
return(1);
}
//.......................... surfaceDraw ............................
// Display rendered surface on display device.
//
short CSurface::surfaceDraw(_HANDLE h, _POINTER iPtr, _POINTER jPtr)
{
// Device dependent. Include code to take lambert shaded bitmap
// and display it on the desired device. For Windows define
// _HANDLE to be CDC*; for X-Windows define _HANDLE as Widget and
// _POINTER as XtPointer.
return(1);
}
//.................................. surfaceFindSmax ................
// Compute intensity scale factor. Smax is found by computing
// the slope of the "best fit" line to the cumulative elevation
// histogram. The equation of an arbitary line is y = 0.5*(1 + x*s)
// where s/2 is the slope. The equation of the fit error is:
//
// E[s] = sum{i: (0.5*(1+i*s) - m_Hz[i])^2 }
// -m_Nz <= i <= m_Nz
//
// Taking the derivative with respect to s gives:
//
// dE[s]/ds = sum {i: i*(0.5*(1+i*s) - m_Hz[i]) }
// -m_Nz <= i <= m_Nz
//
// Setting to zero and solving for s/2 = Smax gives:
//
//Smax = sumH/sumI
//
// where
//
// sumH = sum{i: i*(m_Hz[i] - m_Hz[-i]) }
// 1 <= i <= m_Nz
//
// sumI = sum{i: 2*i*i }
// 1 <= i <= m_Nz
//
short CSurface::surfaceFindSmax(void)
{
double sumI,sumH; // Least squares sums
double hmax,hp,hn; // Histogram values
long i; // Counters
if(surfaceHistogram() == 0) // compute histogram
return(0); // return 0 if unsucessful
hmax = m_Hz[2*m_Nz]; // max elevation difference
sumI = 0.0; // initialize sums
sumH = 0.0;
for(i = 1; i <= m_Nz; i++) // Go through histogram ...
{ // Normalize histogram
hp = ((double) m_Hz[m_Nz+i])/hmax; // values at -/+ i
hn = ((double) m_Hz[m_Nz-i])/hmax; // Compute sums need for
sumH += i*(hp - hn); // Smax calculation
sumI += 2.0*i*i;
}
m_Smax = sumH/sumI; // compute Smax from slope
return(1); // return sucess.
}
//.................................. surfaceHistogram ...............
// Compute elevation difference histogram. A histogram of elevation
// differences in both directions is accumulated for the surface.
// The elevation differences Dx(x,y) and Dy(z,y) are determined form
// adjacent points:
//
// Dx(x,y) = m_Elevations[x+1][y] - m_Elevations[x-1][y]
//
// Dy(x,y) = m_Elevations[x][y+1] - m_Elevations[x][y+1]
//
// Both Dx and Dy are limited to lie between 0 and 2*m_Nz. If
// possible, it is more efficient to compute m_Hz when reading in the
// surface data.
//
short CSurface::surfaceHistogram(void)
{
long i,j; // Counters
long Dx,Dy; // Elevation differences
data _HUGE *Zleft; // Elevation pointers
data _HUGE *Zright;
data _HUGE *Ztop;
data _HUGE *Zbottom;
m_Nz = (long) (m_GridSize + 0.5); // Histogram number
m_Hz = new data[2*m_Nz]; // allocate memory
if(m_Hz == 0) return(0); // no memory for
// computing Smax
for(i = 0; i <= 2*m_Nz; i++) // zero histogram
m_Hz[i] = 0;
Zleft = (m_Elevations + m_Nx + 1); // left elevation
Zright = (m_Elevations + m_Nx + 3); // right elevation
Zbottom = (m_Elevations + 1); // bottom elevation
Ztop = (m_Elevations + 2*(m_Nx + 1)+1); // top elevation
for (j = 1; j < m_Ny; j++) // compute histogram
{
for (i = 1; i < m_Nx; i++)
{
Dx = m_Nz + (long) (*Zright-*Zleft); // x-difference
Dy = m_Nz + (long) (*Ztop-*Zbottom); // y-difference
Dx = (Dx < 0 ? 0 : Dx); // check greater
Dy = (Dy < 0 ? 0 : Dy); // than zero
Dx = (Dx > 2*m_Nz ? 2*m_Nz : Dx); // check less than
Dy = (Dy > 2*m_Nz ? 2*m_Nz : Dy); // upper bound
m_Hz[Dx] += 1; // increment histogram
m_Hz[Dy] += 1; // at these differences
++Zleft; // increment pointers
++Zright;
++Zbottom;
++Ztop;
}
Zleft += 2; // increment pointers
Zright += 2;
Zbottom += 2;
Ztop += 2;
}
for(i = 1; i <= 2*m_Nz; i++) // sum histogram to make
m_Hz[i] = m_Hz[i] + m_Hz[i - 1]; // it cumulative.
return(1); // return success
}
//................................... surfaceLSRender ...............
// Render surface into an image using Lambert Shading. Each surface
// elevation is map directly to a pixel. The intensity of each pixel
// is computed using the Lambert intensity equation:
//
// I(x,y) = Ix[z(x+1,y) - z(x-1,y)] + Iy[z(x,y+1) - z(x, y-1)] + Io
//
// where z = m_Elevations. In order to avoid non-integer arithmetic
// when the surface elevations are integers and scaling factor, FAC,
// is computed and Ix, Iy, Io are scaled accordingly. This results
// in a substantial increase in speed. Pixels at the bitmap edge are
// not computed since elevations in one direction are not available;
// rather these pixel are assigned the value computed for their
// inside neighbor. The result cannot be noticed.
//
short CSurface::surfaceLSRender(void)
{
byte _HUGE *iPtr; // Image pointer
double x,y,z; // Gradients
double sr2; // square root of 2
long Ix, Iy, Io; // Intensity coeffs
long FAC; // Scale factor
long Dx, Dy; // Terrain slopes
long i, j; // Counters
short I; // Pixel intensity
data _HUGE *Zbottom; // Elevation pointers
data _HUGE *Zleft;
data _HUGE *Zright;
data _HUGE *Ztop;
if(surfaceFindSmax() == 0) return(0); // compute Smax
// determine illumination vector
sr2 = sqrt(2.0);
x = m_Smax*sr2*cos(RAD*ILLUM_LON)*cos(RAD*ILLUM_LAT); // Ix
y = m_Smax*sr2*sin(RAD*ILLUM_LON)*cos(RAD*ILLUM_LAT); // Iy
z = ( x < y ? x : y); // z is min illumination
FAC = 1; // intensity scale factor
while( z < 1.0) // choose scale factor so
{ // z is greater than one.
z = 10.0*z; // This allows integer
FAC = 10*FAC; // arithmetic for
} // determining intensities
// mid intensity
Io = (FAC*IMAX_INTENSITY)*sr2*sin(RAD*ILLUM_LAT);/2;
Ix = (long) ((FAC*IMAX_INTENSITY-Io)*x); // scaled sun vector
Iy = (long) ((FAC*IMAX_INTENSITY-Io)*y); // components
// Intensity generation
Zleft = (m_Elevations + m_Nx + 1); // left elevation
Zright = (m_Elevations + m_Nx + 3); // right elevation
Zbottom = (m_Elevations + 1); // bottom elevation
Ztop = (m_Elevations + 2*(m_Nx+1)+1); // top elevation
iPtr = (m_Image + m_Nx + m_Nb + 1); // pointer to image bitmap
for (j = 1; j < m_Ny; j++) // compute intensity...
{
for (i = 1; i < m_Nx; i++)
{
Dx = *Zright - *Zleft; // x-gradient
Dy = *Ztop - *Zbottom; // y-gradient
I = (short) ((Io+Ix*Dx+Iy*Dy)/FAC); // raw intensity
I = (I > 0 ? I : 0); // lower limit
I = (I < IMAX_INTENSITY ? I : IMAX_INTENSITY);
*iPtr = IMAX_INTENSITY - I; // translate
++Zleft; // increment elevation
++Zright; // pointers
++Zbottom;
++Ztop;
++iPtr; // and bitmap pointer
}
iPtr += (m_Nb + 1); // inc bitmap pointer
Zleft += 2; // increment elevation
Zright += 2; // pointers
Zbottom += 2;
Ztop += 2;
}
// At this point the surface image is complete except for the pixel
// that run along its edge. These pixels are filled using adjacent
// values computed above.
for (j = 0; j < m_Ny; j++) // fill first and last
{ // column
*(m_Image + (m_Nx + m_Nb)*j) =
*(m_Image + (m_Nx + m_Nb)*j + 1);
*(m_Image+(m_Nx+m_Nb)*j+m_Nx) =
*(m_Image+(m_Nx+m_Nb)*j+m_Nx-1);
}
for (i = 0; i <= m_Nx; i++) // fill first and last row
{
*(m_Image + i) = *(m_Image+(m_Nx+m_Nb) +i);
*(m_Image+(m_Nx+m_Nb)*m_Ny+i) =
*(m_Image+(m_Nx+m_Nb)*(m_Ny-1)+ i);
}
return(1); // return success
}
//End of File