/src/gdal/frmts/northwood/northwood.cpp

Source (jump to first uncovered line)
/******************************************************************************
 *
 * Project:  GRC/GRD Reader
 * Purpose:  Northwood Format basic implementation
 * Author:   Perry Casson
 *
 ******************************************************************************
 * Copyright (c) 2007, Waypoint Information Technology
 * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "gdal_pam.h"

#include "northwood.h"

#include <algorithm>
#include <cassert>
#include <limits>
#include <string>

int nwt_ParseHeader(NWT_GRID *pGrd, const unsigned char *nwtHeader)
{
    /* double dfTmp; */

    if (nwtHeader[4] == '1')
        pGrd->cFormat = 0x00;  // grd - surface type
    else if (nwtHeader[4] == '8')
        pGrd->cFormat = 0x80;  //  grc classified type

    pGrd->stClassDict = nullptr;

    memcpy(&pGrd->fVersion, &nwtHeader[5], sizeof(pGrd->fVersion));
    CPL_LSBPTR32(&pGrd->fVersion);

    unsigned short usTmp;
    memcpy(&usTmp, &nwtHeader[9], 2);
    CPL_LSBPTR16(&usTmp);
    pGrd->nXSide = static_cast<unsigned int>(usTmp);
    if (pGrd->nXSide == 0)
    {
        memcpy(&pGrd->nXSide, &nwtHeader[128], sizeof(pGrd->nXSide));
        CPL_LSBPTR32(&pGrd->nXSide);
    }
    if (pGrd->nXSide <= 1)
        return FALSE;

    memcpy(&usTmp, &nwtHeader[11], 2);
    CPL_LSBPTR16(&usTmp);
    pGrd->nYSide = static_cast<unsigned int>(usTmp);
    if (pGrd->nYSide == 0)
    {
        memcpy(&pGrd->nYSide, &nwtHeader[132], sizeof(pGrd->nYSide));
        CPL_LSBPTR32(&pGrd->nYSide);
    }

    memcpy(&pGrd->dfMinX, &nwtHeader[13], sizeof(pGrd->dfMinX));
    CPL_LSBPTR64(&pGrd->dfMinX);
    memcpy(&pGrd->dfMaxX, &nwtHeader[21], sizeof(pGrd->dfMaxX));
    CPL_LSBPTR64(&pGrd->dfMaxX);
    memcpy(&pGrd->dfMinY, &nwtHeader[29], sizeof(pGrd->dfMinY));
    CPL_LSBPTR64(&pGrd->dfMinY);
    memcpy(&pGrd->dfMaxY, &nwtHeader[37], sizeof(pGrd->dfMaxY));
    CPL_LSBPTR64(&pGrd->dfMaxY);

    assert(pGrd->nXSide > 1);
    pGrd->dfStepSize = (pGrd->dfMaxX - pGrd->dfMinX) / (pGrd->nXSide - 1);
    /* dfTmp = (pGrd->dfMaxY - pGrd->dfMinY) / (pGrd->nYSide - 1); */

    memcpy(&pGrd->fZMin, &nwtHeader[45], sizeof(pGrd->fZMin));
    CPL_LSBPTR32(&pGrd->fZMin);
    memcpy(&pGrd->fZMax, &nwtHeader[49], sizeof(pGrd->fZMax));
    CPL_LSBPTR32(&pGrd->fZMax);
    memcpy(&pGrd->fZMinScale, &nwtHeader[53], sizeof(pGrd->fZMinScale));
    CPL_LSBPTR32(&pGrd->fZMinScale);
    memcpy(&pGrd->fZMaxScale, &nwtHeader[57], sizeof(pGrd->fZMaxScale));
    CPL_LSBPTR32(&pGrd->fZMaxScale);

    memcpy(&pGrd->cDescription, &nwtHeader[61], sizeof(pGrd->cDescription));
    memcpy(&pGrd->cZUnits, &nwtHeader[93], sizeof(pGrd->cZUnits));

    int i;
    memcpy(&i, &nwtHeader[136], 4);
    CPL_LSBPTR32(&i);

    if (i == 1129336130)
    {  // BMPC
        if (nwtHeader[140] & 0x01)
        {
            pGrd->cHillShadeBrightness = nwtHeader[144];
            pGrd->cHillShadeContrast = nwtHeader[145];
        }
    }

    memcpy(&pGrd->cMICoordSys, &nwtHeader[256], sizeof(pGrd->cMICoordSys));
    pGrd->cMICoordSys[sizeof(pGrd->cMICoordSys) - 1] = '\0';

    pGrd->iZUnits = nwtHeader[512];

    if (nwtHeader[513] & 0x80)
        pGrd->bShowGradient = true;

    if (nwtHeader[513] & 0x40)
        pGrd->bShowHillShade = true;

    if (nwtHeader[513] & 0x20)
        pGrd->bHillShadeExists = true;

    memcpy(&pGrd->iNumColorInflections, &nwtHeader[516], 2);
    CPL_LSBPTR16(&pGrd->iNumColorInflections);

    if (pGrd->iNumColorInflections > 32)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Corrupt header");
        pGrd->iNumColorInflections = 0;
        return FALSE;
    }

    for (i = 0; i < pGrd->iNumColorInflections; i++)
    {
        memcpy(&pGrd->stInflection[i].zVal, &nwtHeader[518 + (7 * i)], 4);
        CPL_LSBPTR32(&pGrd->stInflection[i].zVal);
        memcpy(&pGrd->stInflection[i].r, &nwtHeader[522 + (7 * i)], 1);
        memcpy(&pGrd->stInflection[i].g, &nwtHeader[523 + (7 * i)], 1);
        memcpy(&pGrd->stInflection[i].b, &nwtHeader[524 + (7 * i)], 1);
    }

    memcpy(&pGrd->fHillShadeAzimuth, &nwtHeader[966],
           sizeof(pGrd->fHillShadeAzimuth));
    CPL_LSBPTR32(&pGrd->fHillShadeAzimuth);
    memcpy(&pGrd->fHillShadeAngle, &nwtHeader[970],
           sizeof(pGrd->fHillShadeAngle));
    CPL_LSBPTR32(&pGrd->fHillShadeAngle);

    pGrd->cFormat += nwtHeader[1023];  // the msb for grd/grc was already set

    // there are more types than this - need to build other types for testing
    if (pGrd->cFormat & 0x80)
    {
        if (nwtHeader[1023] == 0)
            pGrd->nBitsPerPixel = 16;
        else
            pGrd->nBitsPerPixel = nwtHeader[1023] * 4;
    }
    else
        pGrd->nBitsPerPixel = nwtHeader[1023] * 8;

    if (pGrd->cFormat & 0x80)  // if is GRC load the Dictionary
    {
        vsi_l_offset nPixels =
            static_cast<vsi_l_offset>(pGrd->nXSide) * pGrd->nYSide;
        unsigned int nBytesPerPixel = pGrd->nBitsPerPixel / 8;
        if (nPixels > 0 &&
            (nBytesPerPixel >
                 std::numeric_limits<vsi_l_offset>::max() / nPixels ||
             nPixels * nBytesPerPixel >
                 std::numeric_limits<vsi_l_offset>::max() - 1024))
        {
            CPLError(CE_Failure, CPLE_FileIO,
                     "Invalid file dimension / bits per pixel");
            return FALSE;
        }
        VSIFSeekL(pGrd->fp, 1024 + nPixels * nBytesPerPixel, SEEK_SET);

        if (!VSIFReadL(&usTmp, 2, 1, pGrd->fp))
        {
            CPLError(CE_Failure, CPLE_FileIO, "Read failure, file short?");
            return FALSE;
        }
        CPL_LSBPTR16(&usTmp);
        pGrd->stClassDict = reinterpret_cast<NWT_CLASSIFIED_DICT *>(
            calloc(1, sizeof(NWT_CLASSIFIED_DICT)));
        if (!pGrd->stClassDict)
        {
            return FALSE;
        }

        pGrd->stClassDict->nNumClassifiedItems = usTmp;

        pGrd->stClassDict->stClassifiedItem =
            reinterpret_cast<NWT_CLASSIFIED_ITEM **>(
                calloc(pGrd->stClassDict->nNumClassifiedItems + 1,
                       sizeof(NWT_CLASSIFIED_ITEM *)));
        if (!pGrd->stClassDict->stClassifiedItem)
        {
            pGrd->stClassDict->nNumClassifiedItems = 0;
            return FALSE;
        }

        // load the dictionary
        for (unsigned int iItem = 0;
             iItem < pGrd->stClassDict->nNumClassifiedItems; iItem++)
        {
            NWT_CLASSIFIED_ITEM *psItem =
                pGrd->stClassDict->stClassifiedItem[iItem] =
                    reinterpret_cast<NWT_CLASSIFIED_ITEM *>(
                        calloc(1, sizeof(NWT_CLASSIFIED_ITEM)));
            if (!psItem)
            {
                return FALSE;
            }

            unsigned char cTmp[256];
            if (!VSIFReadL(&cTmp, 9, 1, pGrd->fp))
            {
                CPLError(CE_Failure, CPLE_FileIO, "Read failure, file short?");
                return FALSE;
            }
            memcpy(&psItem->usPixVal, &cTmp[0], 2);
            CPL_LSBPTR16(&psItem->usPixVal);
            memcpy(&psItem->res1, &cTmp[2], 1);
            memcpy(&psItem->r, &cTmp[3], 1);
            memcpy(&psItem->g, &cTmp[4], 1);
            memcpy(&psItem->b, &cTmp[5], 1);
            memcpy(&psItem->res2, &cTmp[6], 1);
            memcpy(&psItem->usLen, &cTmp[7], 2);
            CPL_LSBPTR16(&psItem->usLen);

            if (psItem->usLen > sizeof(psItem->szClassName) - 1)
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Unexpected long class name, %d characters long - "
                         "unable to read file.",
                         psItem->usLen);
                return FALSE;
            }

            // 0-len class names are possible
            psItem->szClassName[0] = '\0';
            if (psItem->usLen > 0 &&
                !VSIFReadL(&psItem->szClassName, psItem->usLen, 1, pGrd->fp))
                return FALSE;
        }
    }

    return TRUE;
}

// Create a color gradient ranging from ZMin to Zmax using the color
// inflections defined in grid
int nwt_LoadColors(NWT_RGB *pMap, int mapSize, NWT_GRID *pGrd)
{
    int i;
    NWT_RGB sColor;
    int nWarkerMark = 0;

    createIP(0, 255, 255, 255, pMap, &nWarkerMark);
    if (pGrd->iNumColorInflections == 0)
        return 0;

    // If Zmin is less than the 1st inflection use the 1st inflections color to
    // the start of the ramp
    if (pGrd->fZMin <= pGrd->stInflection[0].zVal)
    {
        createIP(1, pGrd->stInflection[0].r, pGrd->stInflection[0].g,
                 pGrd->stInflection[0].b, pMap, &nWarkerMark);
    }
    // find what inflections zmin is between
    for (i = 1; i < pGrd->iNumColorInflections; i++)
    {
        if (pGrd->fZMin < pGrd->stInflection[i].zVal)
        {
            // then we must be between i and i-1
            linearColor(&sColor, &pGrd->stInflection[i - 1],
                        &pGrd->stInflection[i], pGrd->fZMin);
            createIP(1, sColor.r, sColor.g, sColor.b, pMap, &nWarkerMark);
            break;
        }
    }
    // the interesting case of zmin beig higher than the max inflection value
    if (i >= pGrd->iNumColorInflections)
    {
        createIP(1, pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
                 pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
                 pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
                 &nWarkerMark);
        createIP(mapSize - 1,
                 pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
                 pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
                 pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
                 &nWarkerMark);
    }
    else
    {
        int index = 0;
        for (; i < pGrd->iNumColorInflections; i++)
        {
            if (pGrd->fZMax < pGrd->stInflection[i].zVal)
            {
                // then we must be between i and i-1
                linearColor(&sColor, &pGrd->stInflection[i - 1],
                            &pGrd->stInflection[i], pGrd->fZMax);
                index = mapSize - 1;
                createIP(index, sColor.r, sColor.g, sColor.b, pMap,
                         &nWarkerMark);
                break;
            }
            // save the inflections between zmin and zmax
            index =
                static_cast<int>(((pGrd->stInflection[i].zVal - pGrd->fZMin) /
                                  (pGrd->fZMax - pGrd->fZMin)) *
                                 mapSize);

            if (index >= mapSize)
                index = mapSize - 1;
            createIP(index, pGrd->stInflection[i].r, pGrd->stInflection[i].g,
                     pGrd->stInflection[i].b, pMap, &nWarkerMark);
        }
        if (index < mapSize - 1)
            createIP(mapSize - 1,
                     pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
                     pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
                     pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
                     &nWarkerMark);
    }
    return 0;
}

// solve for a color between pIPLow and pIPHigh
void linearColor(NWT_RGB *pRGB, NWT_INFLECTION *pIPLow, NWT_INFLECTION *pIPHigh,
                 float fMid)
{
    if (fMid < pIPLow->zVal)
    {
        pRGB->r = pIPLow->r;
        pRGB->g = pIPLow->g;
        pRGB->b = pIPLow->b;
    }
    else if (fMid > pIPHigh->zVal)
    {
        pRGB->r = pIPHigh->r;
        pRGB->g = pIPHigh->g;
        pRGB->b = pIPHigh->b;
    }
    else
    {
        float scale = (fMid - pIPLow->zVal) / (pIPHigh->zVal - pIPLow->zVal);
        pRGB->r = static_cast<unsigned char>(scale * (pIPHigh->r - pIPLow->r) +
                                             pIPLow->r + 0.5);
        pRGB->g = static_cast<unsigned char>(scale * (pIPHigh->g - pIPLow->g) +
                                             pIPLow->g + 0.5);
        pRGB->b = static_cast<unsigned char>(scale * (pIPHigh->b - pIPLow->b) +
                                             pIPLow->b + 0.5);
    }
}

// insert IP's into the map filling as we go
void createIP(int index, unsigned char r, unsigned char g, unsigned char b,
              NWT_RGB *map, int *pnWarkerMark)
{
    int i;

    if (index == 0)
    {
        map[0].r = r;
        map[0].g = g;
        map[0].b = b;
        *pnWarkerMark = 0;
        return;
    }

    if (index <= *pnWarkerMark)
        return;

    int wm = *pnWarkerMark;

    float rslope =
        static_cast<float>(r - map[wm].r) / static_cast<float>(index - wm);
    float gslope =
        static_cast<float>(g - map[wm].g) / static_cast<float>(index - wm);
    float bslope =
        static_cast<float>(b - map[wm].b) / static_cast<float>(index - wm);
    for (i = wm + 1; i < index; i++)
    {
        map[i].r =
            static_cast<unsigned char>(map[wm].r + ((i - wm) * rslope) + 0.5);
        map[i].g =
            static_cast<unsigned char>(map[wm].g + ((i - wm) * gslope) + 0.5);
        map[i].b =
            static_cast<unsigned char>(map[wm].b + ((i - wm) * bslope) + 0.5);
    }
    map[index].r = r;
    map[index].g = g;
    map[index].b = b;
    *pnWarkerMark = index;
    return;
}

void nwt_HillShade(unsigned char *r, unsigned char *g, unsigned char *b,
                   unsigned char *h)
{
    HLS hls;
    NWT_RGB rgb;
    rgb.r = *r;
    rgb.g = *g;
    rgb.b = *b;
    hls = RGBtoHLS(rgb);
    hls.l = static_cast<short>(hls.l + (*h) * HLSMAX / 256);
    rgb = HLStoRGB(hls);

    *r = rgb.r;
    *g = rgb.g;
    *b = rgb.b;
    return;
}

NWT_GRID *nwtOpenGrid(char *filename)
{
    unsigned char nwtHeader[1024];
    VSILFILE *fp = VSIFOpenL(filename, "rb");

    if (fp == nullptr)
    {
        CPLError(CE_Failure, CPLE_OpenFailed, "Can't open %s", filename);
        return nullptr;
    }

    if (!VSIFReadL(nwtHeader, 1024, 1, fp))
    {
        VSIFCloseL(fp);
        return nullptr;
    }

    if (nwtHeader[0] != 'H' || nwtHeader[1] != 'G' || nwtHeader[2] != 'P' ||
        nwtHeader[3] != 'C')
    {
        VSIFCloseL(fp);
        return nullptr;
    }

    NWT_GRID *pGrd = reinterpret_cast<NWT_GRID *>(calloc(1, sizeof(NWT_GRID)));
    if (!pGrd)
    {
        VSIFCloseL(fp);
        return nullptr;
    }

    if (nwtHeader[4] == '1')
        pGrd->cFormat = 0x00;  // grd - surface type
    else if (nwtHeader[4] == '8')
        pGrd->cFormat = 0x80;  //  grc classified type
    else
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "Unhandled Northwood format type = %0xd", nwtHeader[4]);
        VSIFCloseL(fp);
        free(pGrd);
        return nullptr;
    }

    strncpy(pGrd->szFileName, filename, sizeof(pGrd->szFileName));
    pGrd->szFileName[sizeof(pGrd->szFileName) - 1] = '\0';
    pGrd->fp = fp;
    nwt_ParseHeader(pGrd, nwtHeader);

    return pGrd;
}

// close the file and free the mem
void nwtCloseGrid(NWT_GRID *pGrd)
{
    if ((pGrd->cFormat & 0x80) &&
        pGrd->stClassDict)  // if is GRC - free the Dictionary
    {
        for (unsigned int i = 0; i < pGrd->stClassDict->nNumClassifiedItems;
             i++)
        {
            free(pGrd->stClassDict->stClassifiedItem[i]);
        }
        free(pGrd->stClassDict->stClassifiedItem);
        free(pGrd->stClassDict);
    }
    if (pGrd->fp)
        VSIFCloseL(pGrd->fp);
    free(pGrd);
    return;
}

void nwtPrintGridHeader(NWT_GRID *pGrd)
{
    if (pGrd->cFormat & 0x80)
    {
        printf("\n%s\n\nGrid type is Classified ", pGrd->szFileName); /*ok*/
        if (pGrd->cFormat == 0x81)
            printf("4 bit (Less than 16 Classes)"); /*ok*/
        else if (pGrd->cFormat == 0x82)
            printf("8 bit (Less than 256 Classes)"); /*ok*/
        else if (pGrd->cFormat == 0x84)
            printf("16 bit (Less than 65536 Classes)"); /*ok*/
        else
        {
            printf("GRC - Unhandled Format or Type %d", pGrd->cFormat); /*ok*/
            return;
        }
    }
    else
    {
        printf("\n%s\n\nGrid type is Numeric ", pGrd->szFileName); /*ok*/
        if (pGrd->cFormat == 0x00)
            printf("16 bit (Standard Precision)"); /*ok*/
        else if (pGrd->cFormat == 0x01)
            printf("32 bit (High Precision)"); /*ok*/
        else
        {
            printf("GRD - Unhandled Format or Type %d", pGrd->cFormat); /*ok*/
            return;
        }
    }
    printf("\nDim (x,y) = (%u,%u)", pGrd->nXSide, pGrd->nYSide);     /*ok*/
    printf("\nStep Size = %f", pGrd->dfStepSize);                    /*ok*/
    printf("\nBounds = (%f,%f) (%f,%f)", pGrd->dfMinX, pGrd->dfMinY, /*ok*/
           pGrd->dfMaxX, pGrd->dfMaxY);
    printf("\nCoordinate System = %s", pGrd->cMICoordSys); /*ok*/

    if (!(pGrd->cFormat & 0x80))  // print the numeric specific stuff
    {
        printf("\nMin Z = %f Max Z = %f Z Units = %d \"%s\"", /*ok*/
               pGrd->fZMin, pGrd->fZMax, pGrd->iZUnits, pGrd->cZUnits);

        printf("\n\nDisplay Mode ="); /*ok*/
        if (pGrd->bShowGradient)
            printf(" Color Gradient"); /*ok*/

        if (pGrd->bShowGradient && pGrd->bShowHillShade)
            printf(" and"); /*ok*/

        if (pGrd->bShowHillShade)
            printf(" Hill Shading"); /*ok*/

        for (int i = 0; i < pGrd->iNumColorInflections; i++)
        {
            printf("\nColor Inflection %d - %f (%d,%d,%d)", i + 1, /*ok*/
                   pGrd->stInflection[i].zVal, pGrd->stInflection[i].r,
                   pGrd->stInflection[i].g, pGrd->stInflection[i].b);
        }

        if (pGrd->bHillShadeExists)
        {
            printf("\n\nHill Shade Azumith = %.1f Inclination = %.1f " /*ok*/
                   "Brightness = %d Contrast = %d",
                   pGrd->fHillShadeAzimuth, pGrd->fHillShadeAngle,
                   pGrd->cHillShadeBrightness, pGrd->cHillShadeContrast);
        }
        else
            printf("\n\nNo Hill Shade Data"); /*ok*/
    }
    else  // print the classified specific stuff
    {
        printf("\nNumber of Classes defined = %u", /*ok*/
               pGrd->stClassDict->nNumClassifiedItems);
        for (int i = 0;
             i < static_cast<int>(pGrd->stClassDict->nNumClassifiedItems); i++)
        {
            printf("\n%s - (%d,%d,%d)  Raw = %d  %d %d", /*ok*/
                   pGrd->stClassDict->stClassifiedItem[i]->szClassName,
                   pGrd->stClassDict->stClassifiedItem[i]->r,
                   pGrd->stClassDict->stClassifiedItem[i]->g,
                   pGrd->stClassDict->stClassifiedItem[i]->b,
                   pGrd->stClassDict->stClassifiedItem[i]->usPixVal,
                   pGrd->stClassDict->stClassifiedItem[i]->res1,
                   pGrd->stClassDict->stClassifiedItem[i]->res2);
        }
    }
}

HLS RGBtoHLS(NWT_RGB rgb)
{
    /* get R, G, and B out of DWORD */
    short R = rgb.r;
    short G = rgb.g;
    short B = rgb.b;

    /* calculate lightness */
    unsigned char cMax =
        static_cast<unsigned char>(std::max(std::max(R, G), B));
    unsigned char cMin =
        static_cast<unsigned char>(std::min(std::min(R, G), B));
    HLS hls;
    hls.l = (((cMax + cMin) * HLSMAX) + RGBMAX) / (2 * RGBMAX);

    short Rdelta, Gdelta, Bdelta; /* intermediate value: % of spread from max */
    if (cMax == cMin)
    {                      /* r=g=b --> achromatic case */
        hls.s = 0;         /* saturation */
        hls.h = UNDEFINED; /* hue */
    }
    else
    { /* chromatic case */
        /* saturation */
        if (hls.l <= (HLSMAX / 2))
            hls.s = (((cMax - cMin) * HLSMAX) + ((cMax + cMin) / 2)) /
                    (cMax + cMin);
        else
            hls.s =
                (((cMax - cMin) * HLSMAX) + ((2 * RGBMAX - cMax - cMin) / 2)) /
                (2 * RGBMAX - cMax - cMin);

        /* hue */
        Rdelta =
            (((cMax - R) * (HLSMAX / 6)) + ((cMax - cMin) / 2)) / (cMax - cMin);
        Gdelta =
            (((cMax - G) * (HLSMAX / 6)) + ((cMax - cMin) / 2)) / (cMax - cMin);
        Bdelta =
            (((cMax - B) * (HLSMAX / 6)) + ((cMax - cMin) / 2)) / (cMax - cMin);

        if (R == cMax)
            hls.h = Bdelta - Gdelta;
        else if (G == cMax)
            hls.h = (HLSMAX / 3) + Rdelta - Bdelta;
        else /* B == cMax */
            hls.h = ((2 * HLSMAX) / 3) + Gdelta - Rdelta;

        if (hls.h < 0)
            hls.h += HLSMAX;
        if (hls.h > HLSMAX)
            hls.h -= HLSMAX;
    }
    return hls;
}

/* utility routine for HLStoRGB */
static short HueToRGB(short n1, short n2, short hue)
{
    /* range check: note values passed add/subtract thirds of range */
    if (hue < 0)
        hue += HLSMAX;

    if (hue > HLSMAX)
        hue -= HLSMAX;

    /* return r,g, or b value from this tridrant */
    if (hue < (HLSMAX / 6))
        return n1 + (((n2 - n1) * hue + (HLSMAX / 12)) / (HLSMAX / 6));
    if (hue < (HLSMAX / 2))
        return n2;
    if (hue < ((HLSMAX * 2) / 3))
        return n1 + (((n2 - n1) * (((HLSMAX * 2) / 3) - hue) + (HLSMAX / 12)) /
                     (HLSMAX / 6));
    else
        return n1;
}

NWT_RGB HLStoRGB(HLS hls)
{
    NWT_RGB rgb;

    if (hls.s == 0)
    { /* achromatic case */
        rgb.r = static_cast<unsigned char>((hls.l * RGBMAX) / HLSMAX);
        rgb.g = rgb.r;
        rgb.b = rgb.r;
        if (hls.h != UNDEFINED)
        {
            /* ERROR */
        }
    }
    else
    { /* chromatic case */
        /* set up magic numbers */
        short Magic1, Magic2; /* calculated magic numbers (really!) */
        if (hls.l <= (HLSMAX / 2))
            Magic2 = (hls.l * (HLSMAX + hls.s) + (HLSMAX / 2)) / HLSMAX;
        else
            Magic2 = hls.l + hls.s - ((hls.l * hls.s) + (HLSMAX / 2)) / HLSMAX;
        Magic1 = 2 * hls.l - Magic2;

        /* get RGB, change units from HLSMAX to RGBMAX */
        rgb.r = static_cast<unsigned char>(
            (HueToRGB(Magic1, Magic2, hls.h + (HLSMAX / 3)) * RGBMAX +
             (HLSMAX / 2)) /
            HLSMAX);
        rgb.g = static_cast<unsigned char>(
            (HueToRGB(Magic1, Magic2, hls.h) * RGBMAX + (HLSMAX / 2)) / HLSMAX);
        rgb.b = static_cast<unsigned char>(
            (HueToRGB(Magic1, Magic2, hls.h - (HLSMAX / 3)) * RGBMAX +
             (HLSMAX / 2)) /
            HLSMAX);
    }

    return rgb;
}

Coverage Report

Created: 2025-06-09 08:44

Line	Count	Source (jump to first uncovered line)
1		/******************************************************************************
2		*
3		* Project: GRC/GRD Reader
4		* Purpose: Northwood Format basic implementation
5		* Author: Perry Casson
6		*
7		******************************************************************************
8		* Copyright (c) 2007, Waypoint Information Technology
9		* Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
10		*
11		* SPDX-License-Identifier: MIT
12		****************************************************************************/
13
14		#include "gdal_pam.h"
15
16		#include "northwood.h"
17
18		#include <algorithm>
19		#include <cassert>
20		#include <limits>
21		#include <string>
22
23		int nwt_ParseHeader(NWT_GRID pGrd, const unsigned char nwtHeader)
24	14	{
25		/* double dfTmp; */
26
27	14	if (nwtHeader[4] == '1')
28	5	pGrd->cFormat = 0x00; // grd - surface type
29	9	else if (nwtHeader[4] == '8')
30	9	pGrd->cFormat = 0x80; // grc classified type
31
32	14	pGrd->stClassDict = nullptr;
33
34	14	memcpy(&pGrd->fVersion, &nwtHeader[5], sizeof(pGrd->fVersion));
35	14	CPL_LSBPTR32(&pGrd->fVersion);
36
37	14	unsigned short usTmp;
38	14	memcpy(&usTmp, &nwtHeader[9], 2);
39	14	CPL_LSBPTR16(&usTmp);
40	14	pGrd->nXSide = static_cast<unsigned int>(usTmp);
41	14	if (pGrd->nXSide == 0)
42	0	{
43	0	memcpy(&pGrd->nXSide, &nwtHeader[128], sizeof(pGrd->nXSide));
44	0	CPL_LSBPTR32(&pGrd->nXSide);
45	0	}
46	14	if (pGrd->nXSide <= 1)
47	0	return FALSE;
48
49	14	memcpy(&usTmp, &nwtHeader[11], 2);
50	14	CPL_LSBPTR16(&usTmp);
51	14	pGrd->nYSide = static_cast<unsigned int>(usTmp);
52	14	if (pGrd->nYSide == 0)
53	3	{
54	3	memcpy(&pGrd->nYSide, &nwtHeader[132], sizeof(pGrd->nYSide));
55	3	CPL_LSBPTR32(&pGrd->nYSide);
56	3	}
57
58	14	memcpy(&pGrd->dfMinX, &nwtHeader[13], sizeof(pGrd->dfMinX));
59	14	CPL_LSBPTR64(&pGrd->dfMinX);
60	14	memcpy(&pGrd->dfMaxX, &nwtHeader[21], sizeof(pGrd->dfMaxX));
61	14	CPL_LSBPTR64(&pGrd->dfMaxX);
62	14	memcpy(&pGrd->dfMinY, &nwtHeader[29], sizeof(pGrd->dfMinY));
63	14	CPL_LSBPTR64(&pGrd->dfMinY);
64	14	memcpy(&pGrd->dfMaxY, &nwtHeader[37], sizeof(pGrd->dfMaxY));
65	14	CPL_LSBPTR64(&pGrd->dfMaxY);
66
67	14	assert(pGrd->nXSide > 1);
68	14	pGrd->dfStepSize = (pGrd->dfMaxX - pGrd->dfMinX) / (pGrd->nXSide - 1);
69		/* dfTmp = (pGrd->dfMaxY - pGrd->dfMinY) / (pGrd->nYSide - 1); */
70
71	14	memcpy(&pGrd->fZMin, &nwtHeader[45], sizeof(pGrd->fZMin));
72	14	CPL_LSBPTR32(&pGrd->fZMin);
73	14	memcpy(&pGrd->fZMax, &nwtHeader[49], sizeof(pGrd->fZMax));
74	14	CPL_LSBPTR32(&pGrd->fZMax);
75	14	memcpy(&pGrd->fZMinScale, &nwtHeader[53], sizeof(pGrd->fZMinScale));
76	14	CPL_LSBPTR32(&pGrd->fZMinScale);
77	14	memcpy(&pGrd->fZMaxScale, &nwtHeader[57], sizeof(pGrd->fZMaxScale));
78	14	CPL_LSBPTR32(&pGrd->fZMaxScale);
79
80	14	memcpy(&pGrd->cDescription, &nwtHeader[61], sizeof(pGrd->cDescription));
81	14	memcpy(&pGrd->cZUnits, &nwtHeader[93], sizeof(pGrd->cZUnits));
82
83	14	int i;
84	14	memcpy(&i, &nwtHeader[136], 4);
85	14	CPL_LSBPTR32(&i);
86
87	14	if (i == 1129336130)
88	0	{ // BMPC
89	0	if (nwtHeader[140] & 0x01)
90	0	{
91	0	pGrd->cHillShadeBrightness = nwtHeader[144];
92	0	pGrd->cHillShadeContrast = nwtHeader[145];
93	0	}
94	0	}
95
96	14	memcpy(&pGrd->cMICoordSys, &nwtHeader[256], sizeof(pGrd->cMICoordSys));
97	14	pGrd->cMICoordSys[sizeof(pGrd->cMICoordSys) - 1] = '\0';
98
99	14	pGrd->iZUnits = nwtHeader[512];
100
101	14	if (nwtHeader[513] & 0x80)
102	4	pGrd->bShowGradient = true;
103
104	14	if (nwtHeader[513] & 0x40)
105	3	pGrd->bShowHillShade = true;
106
107	14	if (nwtHeader[513] & 0x20)
108	3	pGrd->bHillShadeExists = true;
109
110	14	memcpy(&pGrd->iNumColorInflections, &nwtHeader[516], 2);
111	14	CPL_LSBPTR16(&pGrd->iNumColorInflections);
112
113	14	if (pGrd->iNumColorInflections > 32)
114	0	{
115	0	CPLError(CE_Failure, CPLE_AppDefined, "Corrupt header");
116	0	pGrd->iNumColorInflections = 0;
117	0	return FALSE;
118	0	}
119
120	23	for (i = 0; i < pGrd->iNumColorInflections; i++)
121	9	{
122	9	memcpy(&pGrd->stInflection[i].zVal, &nwtHeader[518 + (7 * i)], 4);
123	9	CPL_LSBPTR32(&pGrd->stInflection[i].zVal);
124	9	memcpy(&pGrd->stInflection[i].r, &nwtHeader[522 + (7 * i)], 1);
125	9	memcpy(&pGrd->stInflection[i].g, &nwtHeader[523 + (7 * i)], 1);
126	9	memcpy(&pGrd->stInflection[i].b, &nwtHeader[524 + (7 * i)], 1);
127	9	}
128
129	14	memcpy(&pGrd->fHillShadeAzimuth, &nwtHeader[966],
130	14	sizeof(pGrd->fHillShadeAzimuth));
131	14	CPL_LSBPTR32(&pGrd->fHillShadeAzimuth);
132	14	memcpy(&pGrd->fHillShadeAngle, &nwtHeader[970],
133	14	sizeof(pGrd->fHillShadeAngle));
134	14	CPL_LSBPTR32(&pGrd->fHillShadeAngle);
135
136	14	pGrd->cFormat += nwtHeader[1023]; // the msb for grd/grc was already set
137
138		// there are more types than this - need to build other types for testing
139	14	if (pGrd->cFormat & 0x80)
140	9	{
141	9	if (nwtHeader[1023] == 0)
142	3	pGrd->nBitsPerPixel = 16;
143	6	else
144	6	pGrd->nBitsPerPixel = nwtHeader[1023] * 4;
145	9	}
146	5	else
147	5	pGrd->nBitsPerPixel = nwtHeader[1023] * 8;
148
149	14	if (pGrd->cFormat & 0x80) // if is GRC load the Dictionary
150	9	{
151	9	vsi_l_offset nPixels =
152	9	static_cast<vsi_l_offset>(pGrd->nXSide) * pGrd->nYSide;
153	9	unsigned int nBytesPerPixel = pGrd->nBitsPerPixel / 8;
154	9	if (nPixels > 0 &&
155	9	(nBytesPerPixel >
156	9	std::numeric_limits<vsi_l_offset>::max() / nPixels \|\|
157	9	nPixels * nBytesPerPixel >
158	9	std::numeric_limits<vsi_l_offset>::max() - 1024))
159	0	{
160	0	CPLError(CE_Failure, CPLE_FileIO,
161	0	"Invalid file dimension / bits per pixel");
162	0	return FALSE;
163	0	}
164	9	VSIFSeekL(pGrd->fp, 1024 + nPixels * nBytesPerPixel, SEEK_SET);
165
166	9	if (!VSIFReadL(&usTmp, 2, 1, pGrd->fp))
167	0	{
168	0	CPLError(CE_Failure, CPLE_FileIO, "Read failure, file short?");
169	0	return FALSE;
170	0	}
171	9	CPL_LSBPTR16(&usTmp);
172	9	pGrd->stClassDict = reinterpret_cast<NWT_CLASSIFIED_DICT *>(
173	9	calloc(1, sizeof(NWT_CLASSIFIED_DICT)));
174	9	if (!pGrd->stClassDict)
175	0	{
176	0	return FALSE;
177	0	}
178
179	9	pGrd->stClassDict->nNumClassifiedItems = usTmp;
180
181	9	pGrd->stClassDict->stClassifiedItem =
182	9	reinterpret_cast<NWT_CLASSIFIED_ITEM **>(
183	9	calloc(pGrd->stClassDict->nNumClassifiedItems + 1,
184	9	sizeof(NWT_CLASSIFIED_ITEM *)));
185	9	if (!pGrd->stClassDict->stClassifiedItem)
186	0	{
187	0	pGrd->stClassDict->nNumClassifiedItems = 0;
188	0	return FALSE;
189	0	}
190
191		// load the dictionary
192	9	for (unsigned int iItem = 0;
193	16	iItem < pGrd->stClassDict->nNumClassifiedItems; iItem++)
194	10	{
195	10	NWT_CLASSIFIED_ITEM *psItem =
196	10	pGrd->stClassDict->stClassifiedItem[iItem] =
197	10	reinterpret_cast<NWT_CLASSIFIED_ITEM *>(
198	10	calloc(1, sizeof(NWT_CLASSIFIED_ITEM)));
199	10	if (!psItem)
200	0	{
201	0	return FALSE;
202	0	}
203
204	10	unsigned char cTmp[256];
205	10	if (!VSIFReadL(&cTmp, 9, 1, pGrd->fp))
206	0	{
207	0	CPLError(CE_Failure, CPLE_FileIO, "Read failure, file short?");
208	0	return FALSE;
209	0	}
210	10	memcpy(&psItem->usPixVal, &cTmp[0], 2);
211	10	CPL_LSBPTR16(&psItem->usPixVal);
212	10	memcpy(&psItem->res1, &cTmp[2], 1);
213	10	memcpy(&psItem->r, &cTmp[3], 1);
214	10	memcpy(&psItem->g, &cTmp[4], 1);
215	10	memcpy(&psItem->b, &cTmp[5], 1);
216	10	memcpy(&psItem->res2, &cTmp[6], 1);
217	10	memcpy(&psItem->usLen, &cTmp[7], 2);
218	10	CPL_LSBPTR16(&psItem->usLen);
219
220	10	if (psItem->usLen > sizeof(psItem->szClassName) - 1)
221	3	{
222	3	CPLError(CE_Failure, CPLE_AppDefined,
223	3	"Unexpected long class name, %d characters long - "
224	3	"unable to read file.",
225	3	psItem->usLen);
226	3	return FALSE;
227	3	}
228
229		// 0-len class names are possible
230	7	psItem->szClassName[0] = '\0';
231	7	if (psItem->usLen > 0 &&
232	7	!VSIFReadL(&psItem->szClassName, psItem->usLen, 1, pGrd->fp))
233	0	return FALSE;
234	7	}
235	9	}
236
237	11	return TRUE;
238	14	}
239
240		// Create a color gradient ranging from ZMin to Zmax using the color
241		// inflections defined in grid
242		int nwt_LoadColors(NWT_RGB pMap, int mapSize, NWT_GRID pGrd)
243	5	{
244	5	int i;
245	5	NWT_RGB sColor;
246	5	int nWarkerMark = 0;
247
248	5	createIP(0, 255, 255, 255, pMap, &nWarkerMark);
249	5	if (pGrd->iNumColorInflections == 0)
250	2	return 0;
251
252		// If Zmin is less than the 1st inflection use the 1st inflections color to
253		// the start of the ramp
254	3	if (pGrd->fZMin <= pGrd->stInflection[0].zVal)
255	1	{
256	1	createIP(1, pGrd->stInflection[0].r, pGrd->stInflection[0].g,
257	1	pGrd->stInflection[0].b, pMap, &nWarkerMark);
258	1	}
259		// find what inflections zmin is between
260	5	for (i = 1; i < pGrd->iNumColorInflections; i++)
261	3	{
262	3	if (pGrd->fZMin < pGrd->stInflection[i].zVal)
263	1	{
264		// then we must be between i and i-1
265	1	linearColor(&sColor, &pGrd->stInflection[i - 1],
266	1	&pGrd->stInflection[i], pGrd->fZMin);
267	1	createIP(1, sColor.r, sColor.g, sColor.b, pMap, &nWarkerMark);
268	1	break;
269	1	}
270	3	}
271		// the interesting case of zmin beig higher than the max inflection value
272	3	if (i >= pGrd->iNumColorInflections)
273	2	{
274	2	createIP(1, pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
275	2	pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
276	2	pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
277	2	&nWarkerMark);
278	2	createIP(mapSize - 1,
279	2	pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
280	2	pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
281	2	pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
282	2	&nWarkerMark);
283	2	}
284	1	else
285	1	{
286	1	int index = 0;
287	5	for (; i < pGrd->iNumColorInflections; i++)
288	4	{
289	4	if (pGrd->fZMax < pGrd->stInflection[i].zVal)
290	0	{
291		// then we must be between i and i-1
292	0	linearColor(&sColor, &pGrd->stInflection[i - 1],
293	0	&pGrd->stInflection[i], pGrd->fZMax);
294	0	index = mapSize - 1;
295	0	createIP(index, sColor.r, sColor.g, sColor.b, pMap,
296	0	&nWarkerMark);
297	0	break;
298	0	}
299		// save the inflections between zmin and zmax
300	4	index =
301	4	static_cast<int>(((pGrd->stInflection[i].zVal - pGrd->fZMin) /
302	4	(pGrd->fZMax - pGrd->fZMin)) *
303	4	mapSize);
304
305	4	if (index >= mapSize)
306	1	index = mapSize - 1;
307	4	createIP(index, pGrd->stInflection[i].r, pGrd->stInflection[i].g,
308	4	pGrd->stInflection[i].b, pMap, &nWarkerMark);
309	4	}
310	1	if (index < mapSize - 1)
311	0	createIP(mapSize - 1,
312	0	pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
313	0	pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
314	0	pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
315	0	&nWarkerMark);
316	1	}
317	3	return 0;
318	5	}
319
320		// solve for a color between pIPLow and pIPHigh
321		void linearColor(NWT_RGB pRGB, NWT_INFLECTION pIPLow, NWT_INFLECTION *pIPHigh,
322		float fMid)
323	1	{
324	1	if (fMid < pIPLow->zVal)
325	0	{
326	0	pRGB->r = pIPLow->r;
327	0	pRGB->g = pIPLow->g;
328	0	pRGB->b = pIPLow->b;
329	0	}
330	1	else if (fMid > pIPHigh->zVal)
331	0	{
332	0	pRGB->r = pIPHigh->r;
333	0	pRGB->g = pIPHigh->g;
334	0	pRGB->b = pIPHigh->b;
335	0	}
336	1	else
337	1	{
338	1	float scale = (fMid - pIPLow->zVal) / (pIPHigh->zVal - pIPLow->zVal);
339	1	pRGB->r = static_cast<unsigned char>(scale * (pIPHigh->r - pIPLow->r) +
340	1	pIPLow->r + 0.5);
341	1	pRGB->g = static_cast<unsigned char>(scale * (pIPHigh->g - pIPLow->g) +
342	1	pIPLow->g + 0.5);
343	1	pRGB->b = static_cast<unsigned char>(scale * (pIPHigh->b - pIPLow->b) +
344	1	pIPLow->b + 0.5);
345	1	}
346	1	}
347
348		// insert IP's into the map filling as we go
349		void createIP(int index, unsigned char r, unsigned char g, unsigned char b,
350		NWT_RGB map, int pnWarkerMark)
351	15	{
352	15	int i;
353
354	15	if (index == 0)
355	5	{
356	5	map[0].r = r;
357	5	map[0].g = g;
358	5	map[0].b = b;
359	5	*pnWarkerMark = 0;
360	5	return;
361	5	}
362
363	10	if (index <= *pnWarkerMark)
364	1	return;
365
366	9	int wm = *pnWarkerMark;
367
368	9	float rslope =
369	9	static_cast<float>(r - map[wm].r) / static_cast<float>(index - wm);
370	9	float gslope =
371	9	static_cast<float>(g - map[wm].g) / static_cast<float>(index - wm);
372	9	float bslope =
373	9	static_cast<float>(b - map[wm].b) / static_cast<float>(index - wm);
374	12.2k	for (i = wm + 1; i < index; i++)
375	12.2k	{
376	12.2k	map[i].r =
377	12.2k	static_cast<unsigned char>(map[wm].r + ((i - wm) * rslope) + 0.5);
378	12.2k	map[i].g =
379	12.2k	static_cast<unsigned char>(map[wm].g + ((i - wm) * gslope) + 0.5);
380	12.2k	map[i].b =
381	12.2k	static_cast<unsigned char>(map[wm].b + ((i - wm) * bslope) + 0.5);
382	12.2k	}
383	9	map[index].r = r;
384	9	map[index].g = g;
385	9	map[index].b = b;
386	9	*pnWarkerMark = index;
387	9	return;
388	10	}
389
390		void nwt_HillShade(unsigned char r, unsigned char g, unsigned char *b,
391		unsigned char *h)
392	0	{
393	0	HLS hls;
394	0	NWT_RGB rgb;
395	0	rgb.r = *r;
396	0	rgb.g = *g;
397	0	rgb.b = *b;
398	0	hls = RGBtoHLS(rgb);
399	0	hls.l = static_cast<short>(hls.l + (h) HLSMAX / 256);
400	0	rgb = HLStoRGB(hls);
401
402	0	*r = rgb.r;
403	0	*g = rgb.g;
404	0	*b = rgb.b;
405	0	return;
406	0	}
407
408		NWT_GRID nwtOpenGrid(char filename)
409	0	{
410	0	unsigned char nwtHeader[1024];
411	0	VSILFILE *fp = VSIFOpenL(filename, "rb");
412
413	0	if (fp == nullptr)
414	0	{
415	0	CPLError(CE_Failure, CPLE_OpenFailed, "Can't open %s", filename);
416	0	return nullptr;
417	0	}
418
419	0	if (!VSIFReadL(nwtHeader, 1024, 1, fp))
420	0	{
421	0	VSIFCloseL(fp);
422	0	return nullptr;
423	0	}
424
425	0	if (nwtHeader[0] != 'H' \|\| nwtHeader[1] != 'G' \|\| nwtHeader[2] != 'P' \|\|
426	0	nwtHeader[3] != 'C')
427	0	{
428	0	VSIFCloseL(fp);
429	0	return nullptr;
430	0	}
431
432	0	NWT_GRID pGrd = reinterpret_cast<NWT_GRID >(calloc(1, sizeof(NWT_GRID)));
433	0	if (!pGrd)
434	0	{
435	0	VSIFCloseL(fp);
436	0	return nullptr;
437	0	}
438
439	0	if (nwtHeader[4] == '1')
440	0	pGrd->cFormat = 0x00; // grd - surface type
441	0	else if (nwtHeader[4] == '8')
442	0	pGrd->cFormat = 0x80; // grc classified type
443	0	else
444	0	{
445	0	CPLError(CE_Failure, CPLE_NotSupported,
446	0	"Unhandled Northwood format type = %0xd", nwtHeader[4]);
447	0	VSIFCloseL(fp);
448	0	free(pGrd);
449	0	return nullptr;
450	0	}
451
452	0	strncpy(pGrd->szFileName, filename, sizeof(pGrd->szFileName));
453	0	pGrd->szFileName[sizeof(pGrd->szFileName) - 1] = '\0';
454	0	pGrd->fp = fp;
455	0	nwt_ParseHeader(pGrd, nwtHeader);
456
457	0	return pGrd;
458	0	}
459
460		// close the file and free the mem
461		void nwtCloseGrid(NWT_GRID *pGrd)
462	14	{
463	14	if ((pGrd->cFormat & 0x80) &&
464	14	pGrd->stClassDict) // if is GRC - free the Dictionary
465	9	{
466	57.4k	for (unsigned int i = 0; i < pGrd->stClassDict->nNumClassifiedItems;
467	57.4k	i++)
468	57.4k	{
469	57.4k	free(pGrd->stClassDict->stClassifiedItem[i]);
470	57.4k	}
471	9	free(pGrd->stClassDict->stClassifiedItem);
472	9	free(pGrd->stClassDict);
473	9	}
474	14	if (pGrd->fp)
475	0	VSIFCloseL(pGrd->fp);
476	14	free(pGrd);
477	14	return;
478	14	}
479
480		void nwtPrintGridHeader(NWT_GRID *pGrd)
481	0	{
482	0	if (pGrd->cFormat & 0x80)
483	0	{
484	0	printf("\n%s\n\nGrid type is Classified ", pGrd->szFileName); /ok/
485	0	if (pGrd->cFormat == 0x81)
486	0	printf("4 bit (Less than 16 Classes)"); /ok/
487	0	else if (pGrd->cFormat == 0x82)
488	0	printf("8 bit (Less than 256 Classes)"); /ok/
489	0	else if (pGrd->cFormat == 0x84)
490	0	printf("16 bit (Less than 65536 Classes)"); /ok/
491	0	else
492	0	{
493	0	printf("GRC - Unhandled Format or Type %d", pGrd->cFormat); /ok/
494	0	return;
495	0	}
496	0	}
497	0	else
498	0	{
499	0	printf("\n%s\n\nGrid type is Numeric ", pGrd->szFileName); /ok/
500	0	if (pGrd->cFormat == 0x00)
501	0	printf("16 bit (Standard Precision)"); /ok/
502	0	else if (pGrd->cFormat == 0x01)
503	0	printf("32 bit (High Precision)"); /ok/
504	0	else
505	0	{
506	0	printf("GRD - Unhandled Format or Type %d", pGrd->cFormat); /ok/
507	0	return;
508	0	}
509	0	}
510	0	printf("\nDim (x,y) = (%u,%u)", pGrd->nXSide, pGrd->nYSide); /ok/
511	0	printf("\nStep Size = %f", pGrd->dfStepSize); /ok/
512	0	printf("\nBounds = (%f,%f) (%f,%f)", pGrd->dfMinX, pGrd->dfMinY, /ok/
513	0	pGrd->dfMaxX, pGrd->dfMaxY);
514	0	printf("\nCoordinate System = %s", pGrd->cMICoordSys); /ok/
515
516	0	if (!(pGrd->cFormat & 0x80)) // print the numeric specific stuff
517	0	{
518	0	printf("\nMin Z = %f Max Z = %f Z Units = %d \"%s\"", /ok/
519	0	pGrd->fZMin, pGrd->fZMax, pGrd->iZUnits, pGrd->cZUnits);
520
521	0	printf("\n\nDisplay Mode ="); /ok/
522	0	if (pGrd->bShowGradient)
523	0	printf(" Color Gradient"); /ok/
524
525	0	if (pGrd->bShowGradient && pGrd->bShowHillShade)
526	0	printf(" and"); /ok/
527
528	0	if (pGrd->bShowHillShade)
529	0	printf(" Hill Shading"); /ok/
530
531	0	for (int i = 0; i < pGrd->iNumColorInflections; i++)
532	0	{
533	0	printf("\nColor Inflection %d - %f (%d,%d,%d)", i + 1, /ok/
534	0	pGrd->stInflection[i].zVal, pGrd->stInflection[i].r,
535	0	pGrd->stInflection[i].g, pGrd->stInflection[i].b);
536	0	}
537
538	0	if (pGrd->bHillShadeExists)
539	0	{
540	0	printf("\n\nHill Shade Azumith = %.1f Inclination = %.1f " /ok/
541	0	"Brightness = %d Contrast = %d",
542	0	pGrd->fHillShadeAzimuth, pGrd->fHillShadeAngle,
543	0	pGrd->cHillShadeBrightness, pGrd->cHillShadeContrast);
544	0	}
545	0	else
546	0	printf("\n\nNo Hill Shade Data"); /ok/
547	0	}
548	0	else // print the classified specific stuff
549	0	{
550	0	printf("\nNumber of Classes defined = %u", /ok/
551	0	pGrd->stClassDict->nNumClassifiedItems);
552	0	for (int i = 0;
553	0	i < static_cast<int>(pGrd->stClassDict->nNumClassifiedItems); i++)
554	0	{
555	0	printf("\n%s - (%d,%d,%d) Raw = %d %d %d", /ok/
556	0	pGrd->stClassDict->stClassifiedItem[i]->szClassName,
557	0	pGrd->stClassDict->stClassifiedItem[i]->r,
558	0	pGrd->stClassDict->stClassifiedItem[i]->g,
559	0	pGrd->stClassDict->stClassifiedItem[i]->b,
560	0	pGrd->stClassDict->stClassifiedItem[i]->usPixVal,
561	0	pGrd->stClassDict->stClassifiedItem[i]->res1,
562	0	pGrd->stClassDict->stClassifiedItem[i]->res2);
563	0	}
564	0	}
565	0	}
566
567		HLS RGBtoHLS(NWT_RGB rgb)
568	0	{
569		/* get R, G, and B out of DWORD */
570	0	short R = rgb.r;
571	0	short G = rgb.g;
572	0	short B = rgb.b;
573
574		/* calculate lightness */
575	0	unsigned char cMax =
576	0	static_cast<unsigned char>(std::max(std::max(R, G), B));
577	0	unsigned char cMin =
578	0	static_cast<unsigned char>(std::min(std::min(R, G), B));
579	0	HLS hls;
580	0	hls.l = (((cMax + cMin) * HLSMAX) + RGBMAX) / (2 * RGBMAX);
581
582	0	short Rdelta, Gdelta, Bdelta; /* intermediate value: % of spread from max */
583	0	if (cMax == cMin)
584	0	{ /* r=g=b --> achromatic case */
585	0	hls.s = 0; /* saturation */
586	0	hls.h = UNDEFINED; /* hue */
587	0	}
588	0	else
589	0	{ /* chromatic case */
590		/* saturation */
591	0	if (hls.l <= (HLSMAX / 2))
592	0	hls.s = (((cMax - cMin) * HLSMAX) + ((cMax + cMin) / 2)) /
593	0	(cMax + cMin);
594	0	else
595	0	hls.s =
596	0	(((cMax - cMin) * HLSMAX) + ((2 * RGBMAX - cMax - cMin) / 2)) /
597	0	(2 * RGBMAX - cMax - cMin);
598
599		/* hue */
600	0	Rdelta =
601	0	(((cMax - R) * (HLSMAX / 6)) + ((cMax - cMin) / 2)) / (cMax - cMin);
602	0	Gdelta =
603	0	(((cMax - G) * (HLSMAX / 6)) + ((cMax - cMin) / 2)) / (cMax - cMin);
604	0	Bdelta =
605	0	(((cMax - B) * (HLSMAX / 6)) + ((cMax - cMin) / 2)) / (cMax - cMin);
606
607	0	if (R == cMax)
608	0	hls.h = Bdelta - Gdelta;
609	0	else if (G == cMax)
610	0	hls.h = (HLSMAX / 3) + Rdelta - Bdelta;
611	0	else /* B == cMax */
612	0	hls.h = ((2 * HLSMAX) / 3) + Gdelta - Rdelta;
613
614	0	if (hls.h < 0)
615	0	hls.h += HLSMAX;
616	0	if (hls.h > HLSMAX)
617	0	hls.h -= HLSMAX;
618	0	}
619	0	return hls;
620	0	}
621
622		/* utility routine for HLStoRGB */
623		static short HueToRGB(short n1, short n2, short hue)
624	0	{
625		/* range check: note values passed add/subtract thirds of range */
626	0	if (hue < 0)
627	0	hue += HLSMAX;
628
629	0	if (hue > HLSMAX)
630	0	hue -= HLSMAX;
631
632		/* return r,g, or b value from this tridrant */
633	0	if (hue < (HLSMAX / 6))
634	0	return n1 + (((n2 - n1) * hue + (HLSMAX / 12)) / (HLSMAX / 6));
635	0	if (hue < (HLSMAX / 2))
636	0	return n2;
637	0	if (hue < ((HLSMAX * 2) / 3))
638	0	return n1 + (((n2 - n1) * (((HLSMAX * 2) / 3) - hue) + (HLSMAX / 12)) /
639	0	(HLSMAX / 6));
640	0	else
641	0	return n1;
642	0	}
643
644		NWT_RGB HLStoRGB(HLS hls)
645	0	{
646	0	NWT_RGB rgb;
647
648	0	if (hls.s == 0)
649	0	{ /* achromatic case */
650	0	rgb.r = static_cast<unsigned char>((hls.l * RGBMAX) / HLSMAX);
651	0	rgb.g = rgb.r;
652	0	rgb.b = rgb.r;
653	0	if (hls.h != UNDEFINED)
654	0	{
655		/* ERROR */
656	0	}
657	0	}
658	0	else
659	0	{ /* chromatic case */
660		/* set up magic numbers */
661	0	short Magic1, Magic2; /* calculated magic numbers (really!) */
662	0	if (hls.l <= (HLSMAX / 2))
663	0	Magic2 = (hls.l * (HLSMAX + hls.s) + (HLSMAX / 2)) / HLSMAX;
664	0	else
665	0	Magic2 = hls.l + hls.s - ((hls.l * hls.s) + (HLSMAX / 2)) / HLSMAX;
666	0	Magic1 = 2 * hls.l - Magic2;
667
668		/* get RGB, change units from HLSMAX to RGBMAX */
669	0	rgb.r = static_cast<unsigned char>(
670	0	(HueToRGB(Magic1, Magic2, hls.h + (HLSMAX / 3)) * RGBMAX +
671	0	(HLSMAX / 2)) /
672	0	HLSMAX);
673	0	rgb.g = static_cast<unsigned char>(
674	0	(HueToRGB(Magic1, Magic2, hls.h) * RGBMAX + (HLSMAX / 2)) / HLSMAX);
675	0	rgb.b = static_cast<unsigned char>(
676	0	(HueToRGB(Magic1, Magic2, hls.h - (HLSMAX / 3)) * RGBMAX +
677	0	(HLSMAX / 2)) /
678	0	HLSMAX);
679	0	}
680
681	0	return rgb;
682	0	}