overview.cpp

Revision 9398, 29.0 kB (checked in by fwarmerdam, 3 years ago)
updated contact info
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`

Line
1	/******************************************************************************
2	* $Id$
3	*
4	* Project: GDAL Core
5	* Purpose: Helper code to implement overview support in different drivers.
6	* Author: Frank Warmerdam, warmerdam@pobox.com
7	*
8	******************************************************************************
9	* Copyright (c) 2000, Frank Warmerdam
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	* $Log$
31	* Revision 1.14 2006/03/28 14:49:56 fwarmerdam
32	* updated contact info
33	*
34	* Revision 1.13 2006/02/07 19:11:36 fwarmerdam
35	* applied some strategic improved outofmemory checking
36	*
37	* Revision 1.12 2005/04/04 15:24:48 fwarmerdam
38	* Most C entry points now CPL_STDCALL
39	*
40	* Revision 1.11 2003/05/21 04:24:38 warmerda
41	* avoid warnings with cast
42	*
43	* Revision 1.10 2002/07/09 20:33:12 warmerda
44	* expand tabs
45	*
46	* Revision 1.9 2002/04/16 17:49:29 warmerda
47	* Ensure dfRatio is assigned a value.
48	*
49	* Revision 1.8 2001/07/18 04:04:30 warmerda
50	* added CPL_CVSID
51	*
52	* Revision 1.7 2001/07/16 15:21:46 warmerda
53	* added AVERAGE_MAGPHASE option for complex images
54	*
55	* Revision 1.6 2001/01/30 22:32:42 warmerda
56	* added AVERAGE_MP (magnitude preserving averaging) overview resampling type
57	*
58	* Revision 1.5 2000/11/22 18:41:45 warmerda
59	* fixed bug in complex overview generation
60	*
61	* Revision 1.4 2000/08/18 15:25:06 warmerda
62	* added cascading overview regeneration to speed up averaged overviews
63	*
64	* Revision 1.3 2000/07/17 17:08:45 warmerda
65	* added support for complex data
66	*
67	* Revision 1.2 2000/06/26 22:17:58 warmerda
68	* added scaled progress support
69	*
70	* Revision 1.1 2000/04/21 21:54:05 warmerda
71	* New
72	*
73	*/
74
75	#include "gdal_priv.h"
76
77	CPL_CVSID("$Id$");
78
79	/**********************************************************************/
80	/ GDALDownsampleChunk32R() /
81	/**********************************************************************/
82
83	static CPLErr
84	GDALDownsampleChunk32R( int nSrcWidth, int nSrcHeight,
85	float * pafChunk, int nChunkYOff, int nChunkYSize,
86	GDALRasterBand * poOverview,
87	const char * pszResampling )
88
89	{
90	int nDstYOff, nDstYOff2, nOXSize, nOYSize;
91	float *pafDstScanline;
92
93	nOXSize = poOverview->GetXSize();
94	nOYSize = poOverview->GetYSize();
95
96	pafDstScanline = (float ) VSIMalloc(nOXSize sizeof(float));
97	if( pafDstScanline == NULL )
98	{
99	CPLError( CE_Failure, CPLE_OutOfMemory,
100	"GDALDownsampleChunk32R: Out of memory for line buffer." );
101	return CE_Failure;
102	}
103
104	/ -------------------------------------------------------------------- /
105	/ Figure out the line to start writing to, and the first line /
106	/ to not write to. In theory this approach should ensure that /
107	/ every output line will be written if all input chunks are /
108	/ processed. /
109	/ -------------------------------------------------------------------- /
110	nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
111	nDstYOff2 = (int)
112	(0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
113
114	if( nChunkYOff + nChunkYSize == nSrcHeight )
115	nDstYOff2 = nOYSize;
116
117	/ ==================================================================== /
118	/ Loop over destination scanlines. /
119	/ ==================================================================== /
120	for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
121	{
122	float *pafSrcScanline;
123	int nSrcYOff, nSrcYOff2, iDstPixel;
124
125	nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
126	if( nSrcYOff < nChunkYOff )
127	nSrcYOff = nChunkYOff;
128
129	nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
130	if( nSrcYOff2 > nSrcHeight \|\| iDstLine == nOYSize-1 )
131	nSrcYOff2 = nSrcHeight;
132	if( nSrcYOff2 > nChunkYOff + nChunkYSize )
133	nSrcYOff2 = nChunkYOff + nChunkYSize;
134
135	pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth);
136
137	/ -------------------------------------------------------------------- /
138	/ Loop over destination pixels /
139	/ -------------------------------------------------------------------- /
140	for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
141	{
142	int nSrcXOff, nSrcXOff2;
143
144	nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
145	nSrcXOff2 = (int)
146	(0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
147	if( nSrcXOff2 > nSrcWidth )
148	nSrcXOff2 = nSrcWidth;
149
150	if( EQUALN(pszResampling,"NEAR",4) )
151	{
152	pafDstScanline[iDstPixel] = pafSrcScanline[nSrcXOff];
153	}
154	else if( EQUALN(pszResampling,"AVER",4) )
155	{
156	double dfTotal = 0.0;
157	int nCount = 0, iX, iY;
158
159	for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
160	{
161	for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
162	{
163	dfTotal += pafSrcScanline[iX+(iY-nSrcYOff)*nSrcWidth];
164	nCount++;
165	}
166	}
167
168	CPLAssert( nCount > 0 );
169	if( nCount == 0 )
170	pafDstScanline[iDstPixel] = 0.0;
171	else
172	pafDstScanline[iDstPixel] = (float) (dfTotal / nCount);
173	}
174	}
175
176	poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
177	pafDstScanline, nOXSize, 1, GDT_Float32,
178	0, 0 );
179	}
180
181	CPLFree( pafDstScanline );
182
183	return CE_None;
184	}
185
186	/**********************************************************************/
187	/ GDALDownsampleChunkC32R() /
188	/**********************************************************************/
189
190	static CPLErr
191	GDALDownsampleChunkC32R( int nSrcWidth, int nSrcHeight,
192	float * pafChunk, int nChunkYOff, int nChunkYSize,
193	GDALRasterBand * poOverview,
194	const char * pszResampling )
195
196	{
197	int nDstYOff, nDstYOff2, nOXSize, nOYSize;
198	float *pafDstScanline;
199
200	nOXSize = poOverview->GetXSize();
201	nOYSize = poOverview->GetYSize();
202
203	pafDstScanline = (float ) VSIMalloc(nOXSize sizeof(float) * 2);
204	if( pafDstScanline == NULL )
205	{
206	CPLError( CE_Failure, CPLE_OutOfMemory,
207	"GDALDownsampleChunkC32R: Out of memory for line buffer." );
208	return CE_Failure;
209	}
210
211	/ -------------------------------------------------------------------- /
212	/ Figure out the line to start writing to, and the first line /
213	/ to not write to. In theory this approach should ensure that /
214	/ every output line will be written if all input chunks are /
215	/ processed. /
216	/ -------------------------------------------------------------------- /
217	nDstYOff = (int) (0.5 + (nChunkYOff/(double)nSrcHeight) * nOYSize);
218	nDstYOff2 = (int)
219	(0.5 + ((nChunkYOff+nChunkYSize)/(double)nSrcHeight) * nOYSize);
220
221	if( nChunkYOff + nChunkYSize == nSrcHeight )
222	nDstYOff2 = nOYSize;
223
224	/ ==================================================================== /
225	/ Loop over destination scanlines. /
226	/ ==================================================================== /
227	for( int iDstLine = nDstYOff; iDstLine < nDstYOff2; iDstLine++ )
228	{
229	float *pafSrcScanline;
230	int nSrcYOff, nSrcYOff2, iDstPixel;
231
232	nSrcYOff = (int) (0.5 + (iDstLine/(double)nOYSize) * nSrcHeight);
233	if( nSrcYOff < nChunkYOff )
234	nSrcYOff = nChunkYOff;
235
236	nSrcYOff2 = (int) (0.5 + ((iDstLine+1)/(double)nOYSize) * nSrcHeight);
237	if( nSrcYOff2 > nSrcHeight \|\| iDstLine == nOYSize-1 )
238	nSrcYOff2 = nSrcHeight;
239	if( nSrcYOff2 > nChunkYOff + nChunkYSize )
240	nSrcYOff2 = nChunkYOff + nChunkYSize;
241
242	pafSrcScanline = pafChunk + ((nSrcYOff-nChunkYOff) * nSrcWidth) * 2;
243
244	/ -------------------------------------------------------------------- /
245	/ Loop over destination pixels /
246	/ -------------------------------------------------------------------- /
247	for( iDstPixel = 0; iDstPixel < nOXSize; iDstPixel++ )
248	{
249	int nSrcXOff, nSrcXOff2;
250
251	nSrcXOff = (int) (0.5 + (iDstPixel/(double)nOXSize) * nSrcWidth);
252	nSrcXOff2 = (int)
253	(0.5 + ((iDstPixel+1)/(double)nOXSize) * nSrcWidth);
254	if( nSrcXOff2 > nSrcWidth )
255	nSrcXOff2 = nSrcWidth;
256
257	if( EQUALN(pszResampling,"NEAR",4) )
258	{
259	pafDstScanline[iDstPixel2] = pafSrcScanline[nSrcXOff2];
260	pafDstScanline[iDstPixel2+1] = pafSrcScanline[nSrcXOff2+1];
261	}
262	else if( EQUAL(pszResampling,"AVERAGE_MAGPHASE") )
263	{
264	double dfTotalR = 0.0, dfTotalI = 0.0, dfTotalM = 0.0;
265	int nCount = 0, iX, iY;
266
267	for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
268	{
269	for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
270	{
271	double dfR, dfI;
272
273	dfR = pafSrcScanline[iX2+(iY-nSrcYOff)nSrcWidth*2];
274	dfI = pafSrcScanline[iX2+(iY-nSrcYOff)nSrcWidth*2+1];
275	dfTotalR += dfR;
276	dfTotalI += dfI;
277	dfTotalM += sqrt( dfRdfR + dfIdfI );
278	nCount++;
279	}
280	}
281
282	CPLAssert( nCount > 0 );
283	if( nCount == 0 )
284	{
285	pafDstScanline[iDstPixel*2] = 0.0;
286	pafDstScanline[iDstPixel*2+1] = 0.0;
287	}
288	else
289	{
290	double dfM, dfDesiredM, dfRatio=1.0;
291
292	pafDstScanline[iDstPixel2 ] = (float*) (dfTotalR/nCount);
293	pafDstScanline[iDstPixel2+1] = (float*) (dfTotalI/nCount);
294
295	dfM = sqrt(pafDstScanline[iDstPixel2 ]pafDstScanline[iDstPixel*2 ]
296	+ pafDstScanline[iDstPixel2+1]pafDstScanline[iDstPixel*2+1]);
297	dfDesiredM = dfTotalM / nCount;
298	if( dfM != 0.0 )
299	dfRatio = dfDesiredM / dfM;
300
301	pafDstScanline[iDstPixel2 ] = (float) dfRatio;
302	pafDstScanline[iDstPixel2+1] = (float) dfRatio;
303	}
304	}
305	else if( EQUALN(pszResampling,"AVER",4) )
306	{
307	double dfTotalR = 0.0, dfTotalI = 0.0;
308	int nCount = 0, iX, iY;
309
310	for( iY = nSrcYOff; iY < nSrcYOff2; iY++ )
311	{
312	for( iX = nSrcXOff; iX < nSrcXOff2; iX++ )
313	{
314	dfTotalR += pafSrcScanline[iX2+(iY-nSrcYOff)nSrcWidth*2];
315	dfTotalI += pafSrcScanline[iX2+(iY-nSrcYOff)nSrcWidth*2+1];
316	nCount++;
317	}
318	}
319
320	CPLAssert( nCount > 0 );
321	if( nCount == 0 )
322	{
323	pafDstScanline[iDstPixel*2] = 0.0;
324	pafDstScanline[iDstPixel*2+1] = 0.0;
325	}
326	else
327	{
328	pafDstScanline[iDstPixel2 ] = (float*) (dfTotalR/nCount);
329	pafDstScanline[iDstPixel2+1] = (float*) (dfTotalI/nCount);
330	}
331	}
332	}
333
334	poOverview->RasterIO( GF_Write, 0, iDstLine, nOXSize, 1,
335	pafDstScanline, nOXSize, 1, GDT_CFloat32,
336	0, 0 );
337	}
338
339	CPLFree( pafDstScanline );
340
341	return CE_None;
342	}
343
344	/**********************************************************************/
345	/ GDALRegenerateCascadingOverviews() /
346	/ /
347	/ Generate a list of overviews in order from largest to /
348	/ smallest, computing each from the next larger. /
349	/**********************************************************************/
350
351	static CPLErr
352	GDALRegenerateCascadingOverviews(
353	GDALRasterBand poSrcBand, int* nOverviews, GDALRasterBand **papoOvrBands,
354	const char * pszResampling,
355	GDALProgressFunc pfnProgress, void * pProgressData )
356
357	{
358	/ -------------------------------------------------------------------- /
359	/ First, we must put the overviews in order from largest to /
360	/ smallest. /
361	/ -------------------------------------------------------------------- /
362	int i, j;
363
364	for( i = 0; i < nOverviews-1; i++ )
365	{
366	for( j = 0; j < nOverviews - i - 1; j++ )
367	{
368
369	if( papoOvrBands[j]->GetXSize()
370	* (float) papoOvrBands[j]->GetYSize() <
371	papoOvrBands[j+1]->GetXSize()
372	* (float) papoOvrBands[j+1]->GetYSize() )
373	{
374	GDALRasterBand * poTempBand;
375
376	poTempBand = papoOvrBands[j];
377	papoOvrBands[j] = papoOvrBands[j+1];
378	papoOvrBands[j+1] = poTempBand;
379	}
380	}
381	}
382
383	/ -------------------------------------------------------------------- /
384	/ Count total pixels so we can prepare appropriate scaled /
385	/ progress functions. /
386	/ -------------------------------------------------------------------- /
387	double dfTotalPixels = 0.0;
388
389	for( i = 0; i < nOverviews; i++ )
390	{
391	dfTotalPixels += papoOvrBands[i]->GetXSize()
392	* (double) papoOvrBands[i]->GetYSize();
393	}
394
395	/ -------------------------------------------------------------------- /
396	/ Generate all the bands. /
397	/ -------------------------------------------------------------------- /
398	double dfPixelsProcessed = 0.0;
399
400	for( i = 0; i < nOverviews; i++ )
401	{
402	void *pScaledProgressData;
403	double dfPixels;
404	GDALRasterBand *poBaseBand;
405	CPLErr eErr;
406
407	if( i == 0 )
408	poBaseBand = poSrcBand;
409	else
410	poBaseBand = papoOvrBands[i-1];
411
412	dfPixels = papoOvrBands[i]->GetXSize()
413	* (double) papoOvrBands[i]->GetYSize();
414
415	pScaledProgressData = GDALCreateScaledProgress(
416	dfPixelsProcessed / dfTotalPixels,
417	(dfPixelsProcessed + dfPixels) / dfTotalPixels,
418	pfnProgress, pProgressData );
419
420	eErr = GDALRegenerateOverviews( poBaseBand, 1, papoOvrBands + i,
421	pszResampling,
422	GDALScaledProgress,
423	pScaledProgressData );
424	GDALDestroyScaledProgress( pScaledProgressData );
425
426	if( eErr != CE_None )
427	return eErr;
428
429	dfPixelsProcessed += dfPixels;
430	}
431
432	return CE_None;
433	}
434
435	/**********************************************************************/
436	/ GDALRegenerateOverviews() /
437	/**********************************************************************/
438	CPLErr
439	GDALRegenerateOverviews( GDALRasterBand *poSrcBand,
440	int nOverviews, GDALRasterBand **papoOvrBands,
441	const char * pszResampling,
442	GDALProgressFunc pfnProgress, void * pProgressData )
443
444	{
445	int nFullResYChunk, nWidth;
446	int nFRXBlockSize, nFRYBlockSize;
447	GDALDataType eType;
448
449	/ -------------------------------------------------------------------- /
450	/ If we are operating on multiple overviews, and using /
451	/ averaging, lets do them in cascading order to reduce the /
452	/ amount of computation. /
453	/ -------------------------------------------------------------------- /
454	if( EQUALN(pszResampling,"AVER",4) && nOverviews > 1 )
455	return GDALRegenerateCascadingOverviews( poSrcBand,
456	nOverviews, papoOvrBands,
457	pszResampling,
458	pfnProgress,
459	pProgressData );
460
461	/ -------------------------------------------------------------------- /
462	/ Setup one horizontal swath to read from the raw buffer. /
463	/ -------------------------------------------------------------------- /
464	float *pafChunk;
465
466	poSrcBand->GetBlockSize( &nFRXBlockSize, &nFRYBlockSize );
467
468	if( nFRYBlockSize < 4 \|\| nFRYBlockSize > 256 )
469	nFullResYChunk = 32;
470	else
471	nFullResYChunk = nFRYBlockSize;
472
473	if( GDALDataTypeIsComplex( poSrcBand->GetRasterDataType() ) )
474	eType = GDT_CFloat32;
475	else
476	eType = GDT_Float32;
477
478	nWidth = poSrcBand->GetXSize();
479	pafChunk = (float *)
480	VSIMalloc((GDALGetDataTypeSize(eType)/8) * nFullResYChunk * nWidth );
481
482	if( pafChunk == NULL )
483	{
484	CPLError( CE_Failure, CPLE_OutOfMemory,
485	"Out of memory in GDALRegenerateOverviews()." );
486
487	return CE_Failure;
488	}
489
490	/ -------------------------------------------------------------------- /
491	/ Loop over image operating on chunks. /
492	/ -------------------------------------------------------------------- /
493	int nChunkYOff = 0;
494
495	for( nChunkYOff = 0;
496	nChunkYOff < poSrcBand->GetYSize();
497	nChunkYOff += nFullResYChunk )
498	{
499	if( !pfnProgress( nChunkYOff / (double) poSrcBand->GetYSize(),
500	NULL, pProgressData ) )
501	{
502	CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
503	return CE_Failure;
504	}
505
506	if( nFullResYChunk + nChunkYOff > poSrcBand->GetYSize() )
507	nFullResYChunk = poSrcBand->GetYSize() - nChunkYOff;
508
509	/ read chunk /
510	poSrcBand->RasterIO( GF_Read, 0, nChunkYOff, nWidth, nFullResYChunk,
511	pafChunk, nWidth, nFullResYChunk, eType,
512	0, 0 );
513
514	for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
515	{
516	if( eType == GDT_Float32 )
517	GDALDownsampleChunk32R(nWidth, poSrcBand->GetYSize(),
518	pafChunk, nChunkYOff, nFullResYChunk,
519	papoOvrBands[iOverview], pszResampling);
520	else
521	GDALDownsampleChunkC32R(nWidth, poSrcBand->GetYSize(),
522	pafChunk, nChunkYOff, nFullResYChunk,
523	papoOvrBands[iOverview], pszResampling);
524	}
525	}
526
527	VSIFree( pafChunk );
528
529	/ -------------------------------------------------------------------- /
530	/ Renormalized overview mean / stddev if needed. /
531	/ -------------------------------------------------------------------- /
532	if( EQUAL(pszResampling,"AVERAGE_MP") )
533	{
534	GDALOverviewMagnitudeCorrection( (GDALRasterBandH) poSrcBand,
535	nOverviews,
536	(GDALRasterBandH *) papoOvrBands,
537	GDALDummyProgress, NULL );
538	}
539
540	/ -------------------------------------------------------------------- /
541	/ It can be important to flush out data to overviews. /
542	/ -------------------------------------------------------------------- /
543	for( int iOverview = 0; iOverview < nOverviews; iOverview++ )
544	papoOvrBands[iOverview]->FlushCache();
545
546	pfnProgress( 1.0, NULL, pProgressData );
547
548	return CE_None;
549	}
550
551	/**********************************************************************/
552	/ GDALComputeBandStats() /
553	/**********************************************************************/
554
555	CPLErr CPL_STDCALL
556	GDALComputeBandStats( GDALRasterBandH hSrcBand,
557	int nSampleStep,
558	double pdfMean, double* *pdfStdDev,
559	GDALProgressFunc pfnProgress,
560	void *pProgressData )
561
562	{
563	GDALRasterBand poSrcBand = (GDALRasterBand ) hSrcBand;
564	int iLine, nWidth, nHeight;
565	GDALDataType eType = poSrcBand->GetRasterDataType();
566	GDALDataType eWrkType;
567	int bComplex;
568	float *pafData;
569	double dfSum=0.0, dfSum2=0.0;
570	int nSamples = 0;
571
572	nWidth = poSrcBand->GetXSize();
573	nHeight = poSrcBand->GetYSize();
574
575	if( nSampleStep >= nHeight )
576	nSampleStep = 1;
577
578	bComplex = GDALDataTypeIsComplex(eType);
579	if( bComplex )
580	{
581	pafData = (float ) VSIMalloc(nWidth 2 * sizeof(float));
582	eWrkType = GDT_CFloat32;
583	}
584	else
585	{
586	pafData = (float ) VSIMalloc(nWidth sizeof(float));
587	eWrkType = GDT_Float32;
588	}
589
590	if( pafData == NULL )
591	{
592	CPLError( CE_Failure, CPLE_OutOfMemory,
593	"GDALComputeBandStats: Out of memory for buffer." );
594	return CE_Failure;
595	}
596
597	/ -------------------------------------------------------------------- /
598	/ Loop over all sample lines. /
599	/ -------------------------------------------------------------------- /
600	for( iLine = 0; iLine < nHeight; iLine += nSampleStep )
601	{
602	int iPixel;
603
604	if( !pfnProgress( iLine / (double) nHeight,
605	NULL, pProgressData ) )
606	{
607	CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
608	CPLFree( pafData );
609	return CE_Failure;
610	}
611
612	poSrcBand->RasterIO( GF_Read, 0, iLine, nWidth, 1,
613	pafData, nWidth, 1, eWrkType,
614	0, 0 );
615
616	for( iPixel = 0; iPixel < nWidth; iPixel++ )
617	{
618	float fValue;
619
620	if( bComplex )
621	{
622	// Compute the magnitude of the complex value.
623
624	fValue = (float)
625	sqrt(pafData[iPixel2 ] pafData[iPixel*2 ]
626	+ pafData[iPixel2+1] pafData[iPixel*2+1]);
627	}
628	else
629	{
630	fValue = pafData[iPixel];
631	}
632
633	dfSum += fValue;
634	dfSum2 += fValue * fValue;
635	}
636
637	nSamples += nWidth;
638	}
639
640	if( !pfnProgress( 1.0, NULL, pProgressData ) )
641	{
642	CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
643	CPLFree( pafData );
644	return CE_Failure;
645	}
646
647	/ -------------------------------------------------------------------- /
648	/ Produce the result values. /
649	/ -------------------------------------------------------------------- /
650	if( pdfMean != NULL )
651	*pdfMean = dfSum / nSamples;
652
653	if( pdfStdDev != NULL )
654	{
655	double dfMean = dfSum / nSamples;
656
657	pdfStdDev = sqrt((dfSum2 / nSamples) - (dfMean dfMean));
658	}
659
660	CPLFree( pafData );
661
662	return CE_None;
663	}
664
665	/**********************************************************************/
666	/ GDALOverviewMagnitudeCorrection() /
667	/ /
668	/ Correct the mean and standard deviation of the overviews of /
669	/ the given band to match the base layer approximately. /
670	/**********************************************************************/
671
672	CPLErr
673	GDALOverviewMagnitudeCorrection( GDALRasterBandH hBaseBand,
674	int nOverviewCount,
675	GDALRasterBandH *pahOverviews,
676	GDALProgressFunc pfnProgress,
677	void *pProgressData )
678
679	{
680	CPLErr eErr;
681	double dfOrigMean, dfOrigStdDev;
682
683	/ -------------------------------------------------------------------- /
684	/ Compute mean/stddev for source raster. /
685	/ -------------------------------------------------------------------- /
686	eErr = GDALComputeBandStats( hBaseBand, 2, &dfOrigMean, &dfOrigStdDev,
687	pfnProgress, pProgressData );
688
689	if( eErr != CE_None )
690	return eErr;
691
692	/ -------------------------------------------------------------------- /
693	/ Loop on overview bands. /
694	/ -------------------------------------------------------------------- /
695	int iOverview;
696
697	for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )
698	{
699	GDALRasterBand poOverview = (GDALRasterBand )pahOverviews[iOverview];
700	double dfOverviewMean, dfOverviewStdDev;
701	double dfGain;
702
703	eErr = GDALComputeBandStats( pahOverviews[iOverview], 1,
704	&dfOverviewMean, &dfOverviewStdDev,
705	pfnProgress, pProgressData );
706
707	if( eErr != CE_None )
708	return eErr;
709
710	if( dfOrigStdDev < 0.0001 )
711	dfGain = 1.0;
712	else
713	dfGain = dfOrigStdDev / dfOverviewStdDev;
714
715	/ -------------------------------------------------------------------- /
716	/ Apply gain and offset. /
717	/ -------------------------------------------------------------------- /
718	GDALDataType eWrkType, eType = poOverview->GetRasterDataType();
719	int iLine, nWidth, nHeight, bComplex;
720	float *pafData;
721
722	nWidth = poOverview->GetXSize();
723	nHeight = poOverview->GetYSize();
724
725	bComplex = GDALDataTypeIsComplex(eType);
726	if( bComplex )
727	{
728	pafData = (float ) CPLMalloc(nWidth 2 * sizeof(float));
729	eWrkType = GDT_CFloat32;
730	}
731	else
732	{
733	pafData = (float ) CPLMalloc(nWidth sizeof(float));
734	eWrkType = GDT_Float32;
735	}
736
737	if( pafData == NULL )
738	{
739	CPLError( CE_Failure, CPLE_OutOfMemory,
740	"GDALOverviewMagnitudeCorrection: Out of memory for buffer." );
741	return CE_Failure;
742	}
743
744	for( iLine = 0; iLine < nHeight; iLine++ )
745	{
746	int iPixel;
747
748	if( !pfnProgress( iLine / (double) nHeight,
749	NULL, pProgressData ) )
750	{
751	CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
752	CPLFree( pafData );
753	return CE_Failure;
754	}
755
756	poOverview->RasterIO( GF_Read, 0, iLine, nWidth, 1,
757	pafData, nWidth, 1, eWrkType,
758	0, 0 );
759
760	for( iPixel = 0; iPixel < nWidth; iPixel++ )
761	{
762	if( bComplex )
763	{
764	pafData[iPixel2] = (float) dfGain;
765	pafData[iPixel2+1] = (float) dfGain;
766	}
767	else
768	{
769	pafData[iPixel] = (float)
770	((pafData[iPixel]-dfOverviewMean)*dfGain + dfOrigMean);
771
772	}
773	}
774
775	poOverview->RasterIO( GF_Write, 0, iLine, nWidth, 1,
776	pafData, nWidth, 1, eWrkType,
777	0, 0 );
778	}
779
780	if( !pfnProgress( 1.0, NULL, pProgressData ) )
781	{
782	CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
783	CPLFree( pafData );
784	return CE_Failure;
785	}
786
787	CPLFree( pafData );
788	}
789
790	return CE_None;
791	}

Note: See TracBrowser for help on using the browser.

root/tags/gdal_1_3_2/gcore/overview.cpp

Download in other formats: