Drawing an Image as a 3-D Surface

• Download source - 47.99 KB

Introduction

Image processing functions usually involve manipulation of the pixels and colors and displaying the image as a picture in a typical 2D fashion. For scientific applications, sometimes more insight may be gleaned from image features by drawing the image as a 3D surface. This article will describe code that was developed for an image processing application to draw an image as a 3D surface using OpenGL in a dialog window with interactive rotation and zoom.

Background

The image processing application is built as an MFC multi-document interface and utilizes the Microsoft ATL CImage class since it already has the built-in capability for opening and saving images in the most popular formats. The standard CImage member function CImage::Load reads the image file into a bottom-up DIB. However, this can be inconvenient for accessing pixels, so the image is reformatted into a top-down DIB and 32 bit (since most display adapters are now true color) to standardize pixel handling. See Convert32Bit code.

//-----------------------------------------------------------------------
// This function converts 8, 24, & 32 bit DIBs to 32 bit top-down.
//-----------------------------------------------------------------------
void CImagrDoc::Convert32Bit()
{
    if (m_image.GetBPP() < 8) return;

    byte *t, r, g, b;
    int *p, *q, *buf;
    unsigned long i, j, nx, ny;
    RGBQUAD *pRGB = new RGBQUAD[256];   // For GetDIBColorTable()

    nx = m_image.GetWidth();
    ny = m_image.GetHeight();
    unsigned long n = nx * ny;   // No. of pixels

    // Allocate n sized buffer for temp storage
    if (!(buf = (int *)malloc(n * sizeof(int)))) {
       fMessageBox("Error - " __FUNCTION__, MB_ICONERROR,
                    "malloc() error for size: %d", n);
       return;
    }

    BeginWaitCursor();
    switch (m_image.GetBPP()) {
       case 8:
           if (!(i = GetDIBColorTable(m_image.GetDC(), 0, 256, pRGB))) {
               fMessageBox("Error - " __FUNCTION__, MB_ICONERROR,
                          "GetDIBColorTable() error");
               m_image.ReleaseDC();
               goto End;
           }
           //ATLTRACE2("*"__FUNCTION__" GetDIBColorTable(): %d\n", i);
           m_image.ReleaseDC();

           for (j = 0, q = buf; j < ny; j++) {
               t = (byte *) m_image.GetPixelAddress(0, j);
               for (i = 0; i < nx; i++, t++, q++) {
                  r = pRGB[*t].rgbRed;
                  g = pRGB[*t].rgbGreen;
                  b = pRGB[*t].rgbBlue;
                  *q = RGB(b, g, r);    // CImage is BGR
               }
           }
           break;
       case 24:
           for (j = 0, q = buf; j < ny; j++) {
               // Addr. next row (avoids 24 bit offset bottom-up calc.)
               t = (byte *) m_image.GetPixelAddress(0, j);
               for (i = 0; i < nx; i++, t++, q++) {
                  b = *t;        // CImage is BGR
                  g = *(++t);
                  r = *(++t);
                  *q = RGB(b, g, r);
               }
           }
           break;
       case 32:   // Just need to make top-down
           for (j = 0, q = buf; j < ny; j++) {
               // Addr. next row (avoids bottom-up calc.)
               p = (int *) m_image.GetPixelAddress(0, j);
               for (i = 0; i < nx; i++, p++, q++) {
                  *q = *p;
               }
           }
           break;
    }

    // Start a new CImage
    m_image.Destroy();
    if (!m_image.Create(nx, -(int)ny, 32, 0)) {
       fMessageBox("Error - " __FUNCTION__, MB_ICONERROR,
 "Failed bitmap .Create()");
       goto End;
    }
    p = (int *) m_image.GetBits();   // Ptr to new bitmap (top-down DIB)
    memcpy_s(p, n * sizeof(int), buf, n * sizeof(int)); // Copy buf to bitmap
    m_image.ptype = cRGB;        // Update pixel type

End:
    EndWaitCursor();
    free(buf);
}

Once in top-down mode, pixel manipulations are made much simpler. CImage::GetBits() can be used to return a pointer to the first pixel (top-left pixel) and then subsequent pixels are accessed by simply incrementing the pointer in a single for loop.

Pierre Alliez’s code was very helpful on how to setup the device context and pixel format for an OpenGL dialog window (see www.codeproject.com/KB/openGL/wrl_viewer.aspx). Four function calls are needed to establish the device context:

1. GetSafeHwnd() to get a window handle
2. GetDC(hWnd) to get the hDC
3. SetWindowPixelFormat(hDC) to setup the pixel format, and
4. wglCreateContext(hDC) to create the OpenGL device context

See the OnInitDialog() and SetWindowPixelFormat() code for details.

Using the Code

The application draws the 3D image in three OpenGL drawing modes by a chosen menu selection:

1. GL_LINE_STRIP (pixel rows drawn as lines)
2. GL_QUADS (4-sided polygons), and
3. GL_TRIANGLES (3-sided polygons)

The image base is drawn in the x-z plane with the pixel values extending as hills or valleys into the y-axis. Pixel positions are specified from the array of row and column indices. Since each pixel has a red, green, and blue (RGB) component, these need to be mapped into a single y-value for drawing. This is done by using a grey scale conversion calculation as shown below (but could also be done by simply calculating the average of the three RGB values).

y = RED(*p)*0.299 +
GRN(*p)*0.587 + BLU(*p)*0.114

RED, GRN, and BLU are macros that separate the color components from the 32 bit integer pixel value that p points to. (Note: Due to the way the CImage class stores the blue and red bits, the standard macros GetRValue and GetBValue do not return the red and blue components, respectively, but instead return the blue and red; that is, RED = GetBValue, BLU = GetRValue, and GRN = GetGValue.)

Before calling glVertex3i, pixel values are converted to the nearest integer by a macro called NINT (negative pixel values are accounted for here for certain imaging applications):

#define NINT(f)
((f >= 0) ? (int)(f + .5) : (int)(f - .5))

High-resolution images, e.g. greater than 1k pixels square, may look overly cluttered with too many lines drawn (depending on the display resolution and window size). The variable m_Res is used to reduce the drawing resolution by simply skipping pixels in columns and rows. Interpolation may be used, but this method is faster and has been found to work acceptably. The variable m_Res is a counter of skipped rows and columns so as it increases the drawing resolution decreases. Menu options handle halving or doubling m_Res (i.e. values typically increase as 1, 2, 4, 8, …), as this has been found to work acceptably in reducing lines quickly.

Drawing is done by calling glVertex3i in a glBegin-glEnd block. The GL_LINE_STRIP mode simply draws each pixel in a row as a connected line segment. GL_QUADS mode draws four pixels (the current pixel (i, j), the one to its right (i+1, j), the one on the row below it (i, j+1), and the one to the right of that one (i+1, j+1) and so calls glVertex3i four times for each loop. (Note: For readability, +1 was used here, but the code actually does + m_Res as explained above). The GL_TRIANGLES mode also requires the same four pixels per loop, but requires six calls to glVertex3i (drawing two sets of three) due to the way the polygons are drawn and connected. See the OnPaint() code.

void CPlotimDialog::OnPaint()
{
    CDialog::OnPaint();   			// Base class

    if (!m_pDoc) {
fMessageBox("Error - " __FUNCTION__, MB_ICONERROR, "Invalid document pointer");
       return;
    }

    HWND hWnd = GetSafeHwnd();
    HDC hDC = ::GetDC(hWnd);
    wglMakeCurrent(hDC, m_hGLContext);  	// Select as the rendering context

    CMyImage *image = &m_pDoc->m_image;
    int *p = (int *) image->GetBits();  	// Ptr to bitmap
    unsigned long nx = m_pDoc->GetImageCSize().cx;
    unsigned long ny = m_pDoc->GetImageCSize().cy;
    int r, r2, r3, r4, *q;
    unsigned long i, j;
    int min, max;
    GLfloat xscale, yscale, zscale;

    m_pDoc->GetImageMinMax(&min, &max);
    xscale = m_Scale/(GLfloat)nx;       	// Size may have changed
    zscale = m_Scale/(GLfloat)ny;
    yscale = m_Scale/((GLfloat)(max - min)*(GLfloat)2.); // Normalize

    if (m_New) {   // Begin with a nice rotated view
       glClearColor(m_Bkcolor, m_Bkcolor, m_Bkcolor, 0.0); // Set def. background
       glRotatef(m_xRotation, 0.0, 1.0, 0.0);  	// About y-axis
       glRotatef(m_yRotation, 1.0, 0.0, 0.0);  	// About x-axis
    }
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glPolygonMode(GL_FRONT_AND_BACK, m_Fill);

    glPushMatrix();
    glEnable(GL_DEPTH_TEST); /* Enable hidden-surface-removal */

    glTranslatef(-((GLfloat)nx/(GLfloat)2.) * xscale,
                -((GLfloat)(max + min)/(GLfloat)2.) * yscale,
                 ((GLfloat)ny/(GLfloat)2.) * zscale);
    glScalef(xscale, yscale, zscale);

    //(Add if needed) BeginWaitCursor();
    switch (m_DrawMode) {
       case GL_LINE_STRIP:
           for (j = 0; j < ny; j += m_Res) {
               p = (int *)image->GetPixelAddress(0, j); // Inc. addr. for m_Res
               glBegin(m_DrawMode);
              for (i = 0; i < nx; i += m_Res, p += m_Res) {
                  //r = NINT((RED(*p) + GRN(*p) + BLU(*p))/3.);
                  //(Not much diff.)
                  r = NINT(RED(*p)*0.299 + GRN(*p)*0.587 + BLU(*p)*0.114);
                  MapColor(p);
                  glVertex3i(i, r, -(signed)j);
               }
               glEnd();
           }
           break;
       case GL_QUADS:
           glBegin(m_DrawMode);
           for (j = 0; j < ny; j += m_Res) {
               p = (int *)image->GetPixelAddress(0, j);
               for (i = 0; i < nx; i += m_Res, p += m_Res) {
                  r = NINT(RED(*p)*0.299 + GRN(*p)*0.587 + BLU(*p)*0.114);
                  if (i + m_Res < nx) {
                      r2 = NINT(RED(*(p+m_Res))*0.299
                             + GRN(*(p+m_Res))*0.587
                             + BLU(*(p+m_Res))*0.114);
                  }
                  else break;

                  if (j + m_Res < ny) {
                      q = p + nx*m_Res;
                      r3 = NINT(RED(*q)*0.299
                             + GRN(*q)*0.587
                             + BLU(*q)*0.114);
                      r4 = NINT(RED(*(q+m_Res))*0.299
                             + GRN(*(q+m_Res))*0.587
                             + BLU(*(q+m_Res))*0.114);
                  }
                  else break;

                  MapColor(p);
                  glVertex3i(i, r, -(signed)j);
                  MapColor(p+m_Res);
                  glVertex3i(i + m_Res, r2, -(signed)j);
                  MapColor(q+m_Res);
                  glVertex3i(i + m_Res, r4, -(signed)(j + m_Res));
                  MapColor(q);
                  glVertex3i(i, r3, -(signed)(j + m_Res));
               }
           }
           glEnd();
           break;
       case GL_TRIANGLES:
           glBegin(m_DrawMode);
           for (j = 0; j < ny; j += m_Res) {
               p = (int *)image->GetPixelAddress(0, j);
               for (i = 0; i < nx; i += m_Res, p += m_Res) {
                  r = NINT(RED(*p)*0.299 + GRN(*p)*0.587 + BLU(*p)*0.114);
                  if (i + m_Res < nx) {
                      r2 = NINT(RED(*(p+m_Res))*0.299
                             + GRN(*(p+m_Res))*0.587
                             + BLU(*(p+m_Res))*0.114);
                  }
                  else break;

                  if (j + m_Res < ny) {
                      q = p + nx*m_Res;
                      r3 = NINT(RED(*q)*0.299
                             + GRN(*q)*0.587
                             + BLU(*q)*0.114);
                      r4 = NINT(RED(*(q+m_Res))*0.299
                             + GRN(*(q+m_Res))*0.587
                             + BLU(*(q+m_Res))*0.114);
                  }
                  else break;

                  MapColor(p);
                  glVertex3i(i, r, -(signed)j);
                  MapColor(p+m_Res);
                  glVertex3i(i + m_Res, r2, -(signed)j);
                   MapColor(q);
                  glVertex3i(i, r3, -(signed)(j + m_Res));

                  MapColor(p+m_Res);
                  glVertex3i(i + m_Res, r2, -(signed)j);
                  MapColor(q);
                  glVertex3i(i, r3, -(signed)(j + m_Res));
                  MapColor(q+m_Res);
                  glVertex3i(i + m_Res, r4, -(signed)(j + m_Res));
               }
           }
           glEnd();
           break;
    }
    glPopMatrix();
    glFlush();
    SwapBuffers(hDC);

    if (m_New) {
       TextOut(0, 0, RGB(200,0,0), "Left click+drag to rotate view, "
                                 "right click+drag to zoom");
       m_New = false;
    }

    //EndWaitCursor();
    ::ReleaseDC(hWnd, hDC);
}

The MapColor function handles changing the current drawing color to the passed pixel’s RGB values and makes the call to glColor3f. Pixels with values closer to zero (valleys) appear as darker colored areas in the image and higher valued pixels (hills) as lighter colored regions. The three values passed to glColor3f must be of type float between [0.0, 1.0], so each byte color component (range [0, 255]) is converted to float and divided by the maximum value 255.

glColor3f(RED(*p)/(GLfloat)255.,GRN(*p)/(GLfloat)255.,BLU(*p)/(GLfloat)255.)

By default, the application draws the image as a wireframe mesh (GL_LINE mode), but a menu option changes the variable m_Fill = GL_FILL. This is passed to glPolygonMode to change the mode to filled (solid) polygons. The polygon’s color is determined by the pixels that describe the polygon or line and OpenGL smoothly varies the colors.

Once the image is drawn, it’s very useful to be able to rotate and zoom to view from different angles or even underneath. A left-click and drag is used to rotate the drawing and a right-click and drag is used to zoom. Handling the mouse events is done in the usual way and mouse down coordinates are saved in member variables to pass to OnMouseMove. During the mouse capture event, the OnMouseMove function handles the x, y rotations by calling glRotatef. Zooming is done by changing the scale variable m_Scale which is used in OnPaint to calculate scales for x, y, and z based on the image size and the maximum and minimum pixel values. See OnMouseMove.

//-----------------------------------------------------------------------
// The framework calls this member function when the mouse cursor moves.
//-----------------------------------------------------------------------
void CPlotimDialog::OnMouseMove(UINT nFlags, CPoint point)
{
    if (m_LeftButtonDown) { 		// Left : x / y rotation
       HWND hWnd = GetSafeHwnd();
       HDC hDC = ::GetDC(hWnd);
       wglMakeCurrent(hDC, m_hGLContext);
       glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
       glLoadIdentity();

       m_xRotation += m_LeftDownPos.x - point.x;
       glRotatef(m_xRotation, 0.0, 1.0, 0.0); // About y-axis

       m_yRotation += m_LeftDownPos.y - point.y;
       glRotatef(m_yRotation, 1.0, 0.0, 0.0);	// About x-axis
       //ATLTRACE2("***"__FUNCTION__" (%f, %f)\n", m_xRotation, m_yRotation);

       m_LeftDownPos = point;   		// Save for next call
       InvalidateRect(NULL, FALSE); 		// (Faster) Invalidate();

       ::ReleaseDC(hWnd, hDC);
    }
    else if (m_RightButtonDown) {    	// Right : z translation (zoom)
       m_Scale += (m_RightDownPos.y - point.y)/(float)60.0; // 60.0 zoom rate
       m_RightDownPos = point;
       InvalidateRect(NULL, FALSE);
    }

    CDialog::OnMouseMove(nFlags, point);
}

Conclusion

In conclusion, OpenGL was used in a multidocument MFC based C++ application to draw 2D imagery as a 3D surface. Portions of the application’s code was shown so that it can be used in new applications. Having a 3D surface drawing of complex data and being able to rotate and zoom the drawing interactively can be very useful tools to have in an image processing application.

History

• 1^st June, 2011: Initial version
• 26^th July, 2011: Updated source code - more complete set for users to build a project with