Line data    Source code

       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          14 : int nwt_ParseHeader(NWT_GRID *pGrd, const unsigned char *nwtHeader)
      24             : {
      25             :     /* double dfTmp; */
      26             : 
      27          14 :     if (nwtHeader[4] == '1')
      28          12 :         pGrd->cFormat = 0x00;  // grd - surface type
      29           2 :     else if (nwtHeader[4] == '8')
      30           2 :         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             :     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             :     {
      43           0 :         memcpy(&pGrd->nXSide, &nwtHeader[128], sizeof(pGrd->nXSide));
      44           0 :         CPL_LSBPTR32(&pGrd->nXSide);
      45             :     }
      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             :     {
      54           0 :         memcpy(&pGrd->nYSide, &nwtHeader[132], sizeof(pGrd->nYSide));
      55           0 :         CPL_LSBPTR32(&pGrd->nYSide);
      56             :     }
      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             :     int i;
      84          14 :     memcpy(&i, &nwtHeader[136], 4);
      85          14 :     CPL_LSBPTR32(&i);
      86             : 
      87          14 :     if (i == 1129336130)
      88             :     {  // BMPC
      89           0 :         if (nwtHeader[140] & 0x01)
      90             :         {
      91           0 :             pGrd->cHillShadeBrightness = nwtHeader[144];
      92           0 :             pGrd->cHillShadeContrast = nwtHeader[145];
      93             :         }
      94             :     }
      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           0 :         pGrd->bShowGradient = true;
     103             : 
     104          14 :     if (nwtHeader[513] & 0x40)
     105           0 :         pGrd->bShowHillShade = true;
     106             : 
     107          14 :     if (nwtHeader[513] & 0x20)
     108           0 :         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             :     {
     115           0 :         CPLError(CE_Failure, CPLE_AppDefined, "Corrupt header");
     116           0 :         pGrd->iNumColorInflections = 0;
     117           0 :         return FALSE;
     118             :     }
     119             : 
     120          68 :     for (i = 0; i < pGrd->iNumColorInflections; i++)
     121             :     {
     122          54 :         memcpy(&pGrd->stInflection[i].zVal, &nwtHeader[518 + (7 * i)], 4);
     123          54 :         CPL_LSBPTR32(&pGrd->stInflection[i].zVal);
     124          54 :         memcpy(&pGrd->stInflection[i].r, &nwtHeader[522 + (7 * i)], 1);
     125          54 :         memcpy(&pGrd->stInflection[i].g, &nwtHeader[523 + (7 * i)], 1);
     126          54 :         memcpy(&pGrd->stInflection[i].b, &nwtHeader[524 + (7 * i)], 1);
     127             :     }
     128             : 
     129          14 :     memcpy(&pGrd->fHillShadeAzimuth, &nwtHeader[966],
     130             :            sizeof(pGrd->fHillShadeAzimuth));
     131          14 :     CPL_LSBPTR32(&pGrd->fHillShadeAzimuth);
     132          14 :     memcpy(&pGrd->fHillShadeAngle, &nwtHeader[970],
     133             :            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             :     {
     141           2 :         if (nwtHeader[1023] == 0)
     142           0 :             pGrd->nBitsPerPixel = 16;
     143             :         else
     144           2 :             pGrd->nBitsPerPixel = nwtHeader[1023] * 4;
     145             :     }
     146             :     else
     147          12 :         pGrd->nBitsPerPixel = nwtHeader[1023] * 8;
     148             : 
     149          14 :     if (pGrd->cFormat & 0x80)  // if is GRC load the Dictionary
     150             :     {
     151           2 :         vsi_l_offset nPixels =
     152           2 :             static_cast<vsi_l_offset>(pGrd->nXSide) * pGrd->nYSide;
     153           2 :         unsigned int nBytesPerPixel = pGrd->nBitsPerPixel / 8;
     154           4 :         if (nPixels > 0 &&
     155           2 :             (nBytesPerPixel >
     156           2 :                  std::numeric_limits<vsi_l_offset>::max() / nPixels ||
     157           2 :              nPixels * nBytesPerPixel >
     158           2 :                  std::numeric_limits<vsi_l_offset>::max() - 1024))
     159             :         {
     160           0 :             CPLError(CE_Failure, CPLE_FileIO,
     161             :                      "Invalid file dimension / bits per pixel");
     162           0 :             return FALSE;
     163             :         }
     164           2 :         VSIFSeekL(pGrd->fp, 1024 + nPixels * nBytesPerPixel, SEEK_SET);
     165             : 
     166           2 :         if (!VSIFReadL(&usTmp, 2, 1, pGrd->fp))
     167             :         {
     168           0 :             CPLError(CE_Failure, CPLE_FileIO, "Read failure, file short?");
     169           0 :             return FALSE;
     170             :         }
     171           2 :         CPL_LSBPTR16(&usTmp);
     172           2 :         pGrd->stClassDict = reinterpret_cast<NWT_CLASSIFIED_DICT *>(
     173           2 :             calloc(1, sizeof(NWT_CLASSIFIED_DICT)));
     174           2 :         if (!pGrd->stClassDict)
     175             :         {
     176           0 :             return FALSE;
     177             :         }
     178             : 
     179           2 :         pGrd->stClassDict->nNumClassifiedItems = usTmp;
     180             : 
     181           2 :         pGrd->stClassDict->stClassifiedItem =
     182             :             reinterpret_cast<NWT_CLASSIFIED_ITEM **>(
     183           2 :                 calloc(pGrd->stClassDict->nNumClassifiedItems + 1,
     184             :                        sizeof(NWT_CLASSIFIED_ITEM *)));
     185           2 :         if (!pGrd->stClassDict->stClassifiedItem)
     186             :         {
     187           0 :             pGrd->stClassDict->nNumClassifiedItems = 0;
     188           0 :             return FALSE;
     189             :         }
     190             : 
     191             :         // load the dictionary
     192           8 :         for (unsigned int iItem = 0;
     193           8 :              iItem < pGrd->stClassDict->nNumClassifiedItems; iItem++)
     194             :         {
     195             :             NWT_CLASSIFIED_ITEM *psItem =
     196           6 :                 pGrd->stClassDict->stClassifiedItem[iItem] =
     197             :                     reinterpret_cast<NWT_CLASSIFIED_ITEM *>(
     198           6 :                         calloc(1, sizeof(NWT_CLASSIFIED_ITEM)));
     199           6 :             if (!psItem)
     200             :             {
     201           0 :                 return FALSE;
     202             :             }
     203             : 
     204             :             unsigned char cTmp[256];
     205           6 :             if (!VSIFReadL(&cTmp, 9, 1, pGrd->fp))
     206             :             {
     207           0 :                 CPLError(CE_Failure, CPLE_FileIO, "Read failure, file short?");
     208           0 :                 return FALSE;
     209             :             }
     210           6 :             memcpy(&psItem->usPixVal, &cTmp[0], 2);
     211           6 :             CPL_LSBPTR16(&psItem->usPixVal);
     212           6 :             memcpy(&psItem->res1, &cTmp[2], 1);
     213           6 :             memcpy(&psItem->r, &cTmp[3], 1);
     214           6 :             memcpy(&psItem->g, &cTmp[4], 1);
     215           6 :             memcpy(&psItem->b, &cTmp[5], 1);
     216           6 :             memcpy(&psItem->res2, &cTmp[6], 1);
     217           6 :             memcpy(&psItem->usLen, &cTmp[7], 2);
     218           6 :             CPL_LSBPTR16(&psItem->usLen);
     219             : 
     220           6 :             if (psItem->usLen > sizeof(psItem->szClassName) - 1)
     221             :             {
     222           0 :                 CPLError(CE_Failure, CPLE_AppDefined,
     223             :                          "Unexpected long class name, %d characters long - "
     224             :                          "unable to read file.",
     225           0 :                          psItem->usLen);
     226           0 :                 return FALSE;
     227             :             }
     228             : 
     229             :             // 0-len class names are possible
     230           6 :             psItem->szClassName[0] = '\0';
     231          12 :             if (psItem->usLen > 0 &&
     232           6 :                 !VSIFReadL(&psItem->szClassName, psItem->usLen, 1, pGrd->fp))
     233           0 :                 return FALSE;
     234             :         }
     235             :     }
     236             : 
     237          14 :     return TRUE;
     238             : }
     239             : 
     240             : // Create a color gradient ranging from ZMin to Zmax using the color
     241             : // inflections defined in grid
     242          12 : int nwt_LoadColors(NWT_RGB *pMap, int mapSize, NWT_GRID *pGrd)
     243             : {
     244             :     int i;
     245             :     NWT_RGB sColor;
     246          12 :     int nWarkerMark = 0;
     247             : 
     248          12 :     createIP(0, 255, 255, 255, pMap, &nWarkerMark);
     249          12 :     if (pGrd->iNumColorInflections == 0)
     250           0 :         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          12 :     if (pGrd->fZMin <= pGrd->stInflection[0].zVal)
     255             :     {
     256          12 :         createIP(1, pGrd->stInflection[0].r, pGrd->stInflection[0].g,
     257          12 :                  pGrd->stInflection[0].b, pMap, &nWarkerMark);
     258             :     }
     259             :     // find what inflections zmin is between
     260          12 :     for (i = 1; i < pGrd->iNumColorInflections; i++)
     261             :     {
     262          12 :         if (pGrd->fZMin < pGrd->stInflection[i].zVal)
     263             :         {
     264             :             // then we must be between i and i-1
     265          12 :             linearColor(&sColor, &pGrd->stInflection[i - 1],
     266             :                         &pGrd->stInflection[i], pGrd->fZMin);
     267          12 :             createIP(1, sColor.r, sColor.g, sColor.b, pMap, &nWarkerMark);
     268          12 :             break;
     269             :         }
     270             :     }
     271             :     // the interesting case of zmin beig higher than the max inflection value
     272          12 :     if (i >= pGrd->iNumColorInflections)
     273             :     {
     274           0 :         createIP(1, pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
     275           0 :                  pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
     276           0 :                  pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
     277             :                  &nWarkerMark);
     278           0 :         createIP(mapSize - 1,
     279           0 :                  pGrd->stInflection[pGrd->iNumColorInflections - 1].r,
     280           0 :                  pGrd->stInflection[pGrd->iNumColorInflections - 1].g,
     281           0 :                  pGrd->stInflection[pGrd->iNumColorInflections - 1].b, pMap,
     282             :                  &nWarkerMark);
     283             :     }
     284             :     else
     285             :     {
     286          12 :         int index = 0;
     287          54 :         for (; i < pGrd->iNumColorInflections; i++)
     288             :         {
     289          42 :             if (pGrd->fZMax < pGrd->stInflection[i].zVal)
     290             :             {
     291             :                 // then we must be between i and i-1
     292           0 :                 linearColor(&sColor, &pGrd->stInflection[i - 1],
     293             :                             &pGrd->stInflection[i], pGrd->fZMax);
     294           0 :                 index = mapSize - 1;
     295           0 :                 createIP(index, sColor.r, sColor.g, sColor.b, pMap,
     296             :                          &nWarkerMark);
     297           0 :                 break;
     298             :             }
     299             :             // save the inflections between zmin and zmax
     300          42 :             index =
     301          42 :                 static_cast<int>(((pGrd->stInflection[i].zVal - pGrd->fZMin) /
     302          42 :                                   (pGrd->fZMax - pGrd->fZMin)) *
     303             :                                  mapSize);
     304             : 
     305          42 :             if (index >= mapSize)
     306          12 :                 index = mapSize - 1;
     307          42 :             createIP(index, pGrd->stInflection[i].r, pGrd->stInflection[i].g,
     308          42 :                      pGrd->stInflection[i].b, pMap, &nWarkerMark);
     309             :         }
     310          12 :         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             :                      &nWarkerMark);
     316             :     }
     317          12 :     return 0;
     318             : }
     319             : 
     320             : // solve for a color between pIPLow and pIPHigh
     321          12 : void linearColor(NWT_RGB *pRGB, NWT_INFLECTION *pIPLow, NWT_INFLECTION *pIPHigh,
     322             :                  float fMid)
     323             : {
     324          12 :     if (fMid < pIPLow->zVal)
     325             :     {
     326           0 :         pRGB->r = pIPLow->r;
     327           0 :         pRGB->g = pIPLow->g;
     328           0 :         pRGB->b = pIPLow->b;
     329             :     }
     330          12 :     else if (fMid > pIPHigh->zVal)
     331             :     {
     332           0 :         pRGB->r = pIPHigh->r;
     333           0 :         pRGB->g = pIPHigh->g;
     334           0 :         pRGB->b = pIPHigh->b;
     335             :     }
     336             :     else
     337             :     {
     338          12 :         float scale = (fMid - pIPLow->zVal) / (pIPHigh->zVal - pIPLow->zVal);
     339          12 :         pRGB->r = static_cast<unsigned char>(scale * (pIPHigh->r - pIPLow->r) +
     340          12 :                                              pIPLow->r + 0.5);
     341          12 :         pRGB->g = static_cast<unsigned char>(scale * (pIPHigh->g - pIPLow->g) +
     342          12 :                                              pIPLow->g + 0.5);
     343          12 :         pRGB->b = static_cast<unsigned char>(scale * (pIPHigh->b - pIPLow->b) +
     344          12 :                                              pIPLow->b + 0.5);
     345             :     }
     346          12 : }
     347             : 
     348             : // insert IP's into the map filling as we go
     349          78 : void createIP(int index, unsigned char r, unsigned char g, unsigned char b,
     350             :               NWT_RGB *map, int *pnWarkerMark)
     351             : {
     352             :     int i;
     353             : 
     354          78 :     if (index == 0)
     355             :     {
     356          12 :         map[0].r = r;
     357          12 :         map[0].g = g;
     358          12 :         map[0].b = b;
     359          12 :         *pnWarkerMark = 0;
     360          12 :         return;
     361             :     }
     362             : 
     363          66 :     if (index <= *pnWarkerMark)
     364          12 :         return;
     365             : 
     366          54 :     int wm = *pnWarkerMark;
     367             : 
     368          54 :     float rslope =
     369          54 :         static_cast<float>(r - map[wm].r) / static_cast<float>(index - wm);
     370          54 :     float gslope =
     371          54 :         static_cast<float>(g - map[wm].g) / static_cast<float>(index - wm);
     372          54 :     float bslope =
     373          54 :         static_cast<float>(b - map[wm].b) / static_cast<float>(index - wm);
     374       49140 :     for (i = wm + 1; i < index; i++)
     375             :     {
     376       49086 :         map[i].r =
     377       49086 :             static_cast<unsigned char>(map[wm].r + ((i - wm) * rslope) + 0.5);
     378       49086 :         map[i].g =
     379       49086 :             static_cast<unsigned char>(map[wm].g + ((i - wm) * gslope) + 0.5);
     380       49086 :         map[i].b =
     381       49086 :             static_cast<unsigned char>(map[wm].b + ((i - wm) * bslope) + 0.5);
     382             :     }
     383          54 :     map[index].r = r;
     384          54 :     map[index].g = g;
     385          54 :     map[index].b = b;
     386          54 :     *pnWarkerMark = index;
     387          54 :     return;
     388             : }
     389             : 
     390           0 : void nwt_HillShade(unsigned char *r, unsigned char *g, unsigned char *b,
     391             :                    unsigned char *h)
     392             : {
     393             :     HLS hls;
     394             :     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             : }
     407             : 
     408           0 : NWT_GRID *nwtOpenGrid(char *filename)
     409             : {
     410             :     unsigned char nwtHeader[1024];
     411           0 :     VSILFILE *fp = VSIFOpenL(filename, "rb");
     412             : 
     413           0 :     if (fp == nullptr)
     414             :     {
     415           0 :         CPLError(CE_Failure, CPLE_OpenFailed, "Can't open %s", filename);
     416           0 :         return nullptr;
     417             :     }
     418             : 
     419           0 :     if (!VSIFReadL(nwtHeader, 1024, 1, fp))
     420             :     {
     421           0 :         VSIFCloseL(fp);
     422           0 :         return nullptr;
     423             :     }
     424             : 
     425           0 :     if (nwtHeader[0] != 'H' || nwtHeader[1] != 'G' || nwtHeader[2] != 'P' ||
     426           0 :         nwtHeader[3] != 'C')
     427             :     {
     428           0 :         VSIFCloseL(fp);
     429           0 :         return nullptr;
     430             :     }
     431             : 
     432           0 :     NWT_GRID *pGrd = reinterpret_cast<NWT_GRID *>(calloc(1, sizeof(NWT_GRID)));
     433           0 :     if (!pGrd)
     434             :     {
     435           0 :         VSIFCloseL(fp);
     436           0 :         return nullptr;
     437             :     }
     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             :     else
     444             :     {
     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             :     }
     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             : }
     459             : 
     460             : // close the file and free the mem
     461          28 : void nwtCloseGrid(NWT_GRID *pGrd)
     462             : {
     463          30 :     if ((pGrd->cFormat & 0x80) &&
     464          28 :         pGrd->stClassDict)  // if is GRC - free the Dictionary
     465             :     {
     466           8 :         for (unsigned int i = 0; i < pGrd->stClassDict->nNumClassifiedItems;
     467             :              i++)
     468             :         {
     469           6 :             free(pGrd->stClassDict->stClassifiedItem[i]);
     470             :         }
     471           2 :         free(pGrd->stClassDict->stClassifiedItem);
     472           2 :         free(pGrd->stClassDict);
     473             :     }
     474          28 :     if (pGrd->fp)
     475           0 :         VSIFCloseL(pGrd->fp);
     476          28 :     free(pGrd);
     477          28 :     return;
     478             : }
     479             : 
     480           0 : void nwtPrintGridHeader(NWT_GRID *pGrd)
     481             : {
     482           0 :     if (pGrd->cFormat & 0x80)
     483             :     {
     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             :         else
     492             :         {
     493           0 :             printf("GRC - Unhandled Format or Type %d", pGrd->cFormat); /*ok*/
     494           0 :             return;
     495             :         }
     496             :     }
     497             :     else
     498             :     {
     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             :         else
     505             :         {
     506           0 :             printf("GRD - Unhandled Format or Type %d", pGrd->cFormat); /*ok*/
     507           0 :             return;
     508             :         }
     509             :     }
     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             :            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             :     {
     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             :         {
     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             :         }
     537             : 
     538           0 :         if (pGrd->bHillShadeExists)
     539             :         {
     540           0 :             printf("\n\nHill Shade Azumith = %.1f Inclination = %.1f " /*ok*/
     541             :                    "Brightness = %d Contrast = %d",
     542           0 :                    pGrd->fHillShadeAzimuth, pGrd->fHillShadeAngle,
     543           0 :                    pGrd->cHillShadeBrightness, pGrd->cHillShadeContrast);
     544             :         }
     545             :         else
     546           0 :             printf("\n\nNo Hill Shade Data"); /*ok*/
     547             :     }
     548             :     else  // print the classified specific stuff
     549             :     {
     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             :         {
     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             :         }
     564             :     }
     565             : }
     566             : 
     567           0 : HLS RGBtoHLS(NWT_RGB rgb)
     568             : {
     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             :     unsigned char cMax =
     576           0 :         static_cast<unsigned char>(std::max(std::max(R, G), B));
     577             :     unsigned char cMin =
     578           0 :         static_cast<unsigned char>(std::min(std::min(R, G), B));
     579             :     HLS hls;
     580           0 :     hls.l = (((cMax + cMin) * HLSMAX) + RGBMAX) / (2 * RGBMAX);
     581             : 
     582             :     short Rdelta, Gdelta, Bdelta; /* intermediate value: % of spread from max */
     583           0 :     if (cMax == cMin)
     584             :     {                      /* r=g=b --> achromatic case */
     585           0 :         hls.s = 0;         /* saturation */
     586           0 :         hls.h = UNDEFINED; /* hue */
     587             :     }
     588             :     else
     589             :     { /* 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             :         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             :         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             :     }
     619           0 :     return hls;
     620             : }
     621             : 
     622             : /* utility routine for HLStoRGB */
     623           0 : static short HueToRGB(short n1, short n2, short hue)
     624             : {
     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             :     else
     641           0 :         return n1;
     642             : }
     643             : 
     644           0 : NWT_RGB HLStoRGB(HLS hls)
     645             : {
     646             :     NWT_RGB rgb;
     647             : 
     648           0 :     if (hls.s == 0)
     649             :     { /* 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             :         {
     655             :             /* ERROR */
     656             :         }
     657             :     }
     658             :     else
     659             :     { /* chromatic case */
     660             :         /* set up magic numbers */
     661             :         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             :         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             :             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             :             HLSMAX);
     679             :     }
     680             : 
     681           0 :     return rgb;
     682             : }