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