62 BYTE CSG_Grid::m_Bitmask[8] = { 0x1, 0x2, 0x4, 0x8, 0x10, 0x20, 0x40, 0x80 };
82 if( !pGrid->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
92 if( pGrid->
Create(File, Type, bCached, bLoadData) )
97 delete(pGrid);
return( NULL );
103 pGrid =
new CSG_Grid(pGrid, Type, bCached);
105 if( !pGrid->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
113 if( !pGrid->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
119 CSG_Grid *pGrid =
new CSG_Grid(Type, NX, NY, Cellsize > 0. ? Cellsize : 1., xMin, yMin, bCached);
121 if( !pGrid->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
208 void CSG_Grid::_On_Construction(
void)
214 m_Cache_Stream = NULL;
216 m_Cache_bSwap =
false;
217 m_Cache_bFlip =
false;
247 #pragma omp parallel for
248 for(
int y=0; y<
Get_NY(); y++)
for(
int x=0; x<
Get_NX(); x++)
281 switch( (m_Type = Type) )
301 if( m_System.
Create(System) ==
false )
313 return( _Memory_Create(bCached) );
327 if( _Load_PGSQL (File, bCached, bLoadData)
328 || _Load_Native (File, bCached, bLoadData)
329 || _Load_Compressed(File, bCached, bLoadData)
330 || _Load_Surfer (File, bCached, bLoadData)
331 || _Load_External (File, bCached, bLoadData) )
408 if( (Scale != m_zScale && Scale != 0.) || Offset != m_zOffset )
445 return( m_Values != NULL ||
is_Cached() );
475 return( m_System == System );
495 if( y >= 0 && y <
Get_NY() )
499 for(
int x=0; x<
Get_NX(); x++)
513 for(
int x=0; x<
Get_NX(); x++)
548 return(
Get_Value(p.
x, p.
y, Value, Resampling, bNoData, bByteWise) );
559 if( bNoData ||
is_InGrid(ix + (
int)(0.5 + dx), iy + (
int)(0.5 + dy)) )
582 inline bool CSG_Grid::_Get_ValAtPos_NearestNeighbour(
double &Value,
int x,
int y,
double dx,
double dy)
const
584 if(
is_InGrid(x = x + (
int)(0.5 + dx), y = y + (
int)(0.5 + dy)) )
600 #define BILINEAR_ADD(ix, iy, d) if( is_InGrid(ix, iy) ) {\
601 n += d; z += d * asDouble(ix, iy);\
605 #define BILINEAR_ADD_BYTE(ix, iy, d) if( is_InGrid(ix, iy) ) {\
606 n += d; sLong v = asInt(ix, iy);\
607 z[0] += d * SG_GET_BYTE_0(v);\
608 z[1] += d * SG_GET_BYTE_1(v);\
609 z[2] += d * SG_GET_BYTE_2(v);\
610 z[3] += d * SG_GET_BYTE_3(v);\
614 inline bool CSG_Grid::_Get_ValAtPos_BiLinear(
double &Value,
int x,
int y,
double dx,
double dy,
bool bByteWise)
const
618 double z = 0., n = 0.;
663 inline double CSG_Grid::_Get_ValAtPos_BiCubicSpline(
double dx,
double dy,
double v_xy[4][4])
const
665 #define BiCubicSpline(d, v) (v[1] + 0.5 * d * (v[2] - v[0] + d * (2 * v[0] - 5 * v[1] + 4 * v[2] - v[3] + d * (3 * (v[1] - v[2]) + v[3] - v[0]))))
669 for(
int ix=0; ix<4; ix++)
678 inline bool CSG_Grid::_Get_ValAtPos_BiCubicSpline(
double &Value,
int x,
int y,
double dx,
double dy,
bool bByteWise)
const
684 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
686 Value = _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy);
693 double v_xy[4][4][4];
695 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
698 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[0]),
699 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[1]),
700 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[2]),
701 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[3])
717 inline double CSG_Grid::_Get_ValAtPos_BSpline(
double dx,
double dy,
double v_xy[4][4])
const
721 for(
int i=0; i<4; i++)
726 if( (d = i - dx + 1.) > 0. ) s += d*d*d;
727 if( (d = i - dx + 0.) > 0. ) s += -4. * d*d*d;
728 if( (d = i - dx - 1.) > 0. ) s += 6. * d*d*d;
729 if( (d = i - dx - 2.) > 0. ) s += -4. * d*d*d;
734 if( (d = i - dy + 1.) > 0. ) s += d*d*d;
735 if( (d = i - dy + 0.) > 0. ) s += -4. * d*d*d;
736 if( (d = i - dy - 1.) > 0. ) s += 6. * d*d*d;
737 if( (d = i - dy - 2.) > 0. ) s += -4. * d*d*d;
744 for(
int iy=0; iy<4; iy++)
746 for(
int ix=0; ix<4; ix++)
748 z += v_xy[ix][iy] * Rx[ix] * Ry[iy];
756 inline bool CSG_Grid::_Get_ValAtPos_BSpline(
double &Value,
int x,
int y,
double dx,
double dy,
bool bByteWise)
const
762 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
764 Value = _Get_ValAtPos_BSpline(dx, dy, v_xy);
771 double v_xy[4][4][4];
773 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
776 _Get_ValAtPos_BSpline(dx, dy, v_xy[0]),
777 _Get_ValAtPos_BSpline(dx, dy, v_xy[1]),
778 _Get_ValAtPos_BSpline(dx, dy, v_xy[2]),
779 _Get_ValAtPos_BSpline(dx, dy, v_xy[3])
795 inline bool CSG_Grid::_Get_ValAtPos_Fill4x4Submatrix(
int x,
int y,
double v_xy[4][4])
const
797 int ix, iy, jx, jy, nNoData = 0;
800 for(iy=0, jy=y-1; iy<4; iy++, jy++)
802 for(ix=0, jx=x-1; ix<4; ix++, jx++)
818 for(
int i=0; nNoData>0 && nNoData<16 && i<16; i++)
822 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
824 t_xy[ix][iy] = v_xy[ix][iy];
827 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
834 for(jy=iy-1; jy<=iy+1; jy++)
for(jx=ix-1; jx<=ix+1; jx++)
838 s +=
asDouble(jx + x - 1, jy + y - 1);
841 else if( jy >= 0 && jy < 4 && jx >= 0 && jx < 4 && !
is_NoData_Value(t_xy[jx][jy]) )
850 v_xy[ix][iy] = s / n;
859 return( nNoData == 0 );
863 inline bool CSG_Grid::_Get_ValAtPos_Fill4x4Submatrix(
int x,
int y,
double v_xy[4][4][4])
const
865 int ix, iy, jx, jy, nNoData = 0;
868 for(iy=0, jy=y-1; iy<4; iy++, jy++)
870 for(ix=0, jx=x-1; ix<4; ix++, jx++)
874 int z =
asInt(jx, jy);
883 v_xy[0][ix][iy] = -1;
891 for(
int i=0; nNoData>0 && nNoData<16 && i<16; i++)
893 double t_xy[4][4][4];
895 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
897 t_xy[0][ix][iy] = v_xy[0][ix][iy];
898 t_xy[1][ix][iy] = v_xy[1][ix][iy];
899 t_xy[2][ix][iy] = v_xy[2][ix][iy];
900 t_xy[3][ix][iy] = v_xy[3][ix][iy];
903 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
905 if( t_xy[0][ix][iy] < 0 )
908 double s[4]; s[0] = s[1] = s[2] = s[3] = 0.;
910 for(jy=iy-1; jy<=iy+1; jy++)
for(jx=ix-1; jx<=ix+1; jx++)
914 int z =
asInt(jx + x - 1, jy + y - 1);
922 else if( jy >= 0 && jy < 4 && jx >= 0 && jx < 4 && !
is_NoData_Value(t_xy[0][jx][jy]) )
924 s[0] += t_xy[0][jx][jy];
925 s[1] += t_xy[1][jx][jy];
926 s[2] += t_xy[2][jx][jy];
927 s[3] += t_xy[3][jx][jy];
934 v_xy[0][ix][iy] = s[0] / n;
935 v_xy[1][ix][iy] = s[1] / n;
936 v_xy[2][ix][iy] = s[2] / n;
937 v_xy[3][ix][iy] = s[3] / n;
946 return( nNoData == 0 );
981 m_Statistics += Scaling ? Offset + Scaling * Value : Value;
993 for(
int x=0; x<
Get_NX(); x++)
995 double Value =
asDouble(x, y,
false);
999 m_Statistics += Scaling ? Offset + Scaling * Value : Value;
1055 if( Quantile <= 0. ) {
return(
Get_Min() ); }
1056 if( Quantile >= 1. ) {
return(
Get_Max() ); }
1058 if( bFromHistogram )
1078 return(
Get_Quantile(0.01 * Percentile, bFromHistogram) );
1095 Update();
return( m_Statistics );
1112 if( xMin > xMax || yMin > yMax )
1117 Statistics.
Create(bHoldValues);
1119 int nx = 1 + (xMax - xMin);
1120 int ny = 1 + (yMax - yMin);
1121 sLong nCells = nx * ny;
1129 for(
double i=0; i<(double)nCells; i+=d)
1131 int y = yMin + (int)i / nx;
1132 int x = xMin + (int)i % nx;
1134 double Value =
asDouble(x, y,
false);
1138 Statistics += Scaling ? Offset + Scaling * Value : Value;
1144 for(
int y=yMin; y<=yMax; y++)
1146 for(
int x=xMin; x<=xMax; x++)
1148 double Value =
asDouble(x, y,
false);
1152 Statistics += Scaling ? Offset + Scaling * Value : Value;
1162 #define SG_GRID_HISTOGRAM_CLASSES_DEFAULT 255
1183 return( m_Histogram );
1201 if( xMin > xMax || yMin > yMax )
1208 int nx = 1 + (xMax - xMin);
1209 int ny = 1 + (yMax - yMin);
1210 sLong nCells = nx * ny;
1218 for(
double i=0; i<(double)nCells; i+=d)
1220 int y = yMin + (int)i / nx;
1221 int x = xMin + (int)i % nx;
1223 double Value =
asDouble(x, y,
false);
1227 Histogram += Scaling ? Offset + Scaling * Value : Value;
1233 for(
int y=yMin; y<=yMax; y++)
1235 for(
int x=xMin; x<=xMax; x++)
1237 double Value =
asDouble(x, y,
false);
1241 Histogram += Scaling ? Offset + Scaling * Value : Value;
1247 return( Histogram.
Update() );
1258 bool CSG_Grid::_Set_Index(
void)
1271 sLong i, j, k, l, ir, n, *istack, jstack, nstack, indxt, itemp, nData;
1281 m_Index[--nData] = i;
1320 for(j=l+1; j<=ir; j++)
1325 for(i=j-1; i>=0; i--)
1332 m_Index[i + 1] = m_Index[i];
1335 m_Index[i + 1] = indxt;
1343 ir = istack[jstack--];
1344 l = istack[jstack--];
1350 #define SORT_SWAP(a,b) {itemp=(a);(a)=(b);(b)=itemp;}
1372 do i++;
while(
asDouble(m_Index[i]) < a);
1373 do j--;
while(
asDouble(m_Index[j]) > a);
1385 m_Index[l] = m_Index[j];
1389 if( jstack >= nstack )
1395 if( ir - i + 1 >= j - l )
1397 istack[jstack] = ir;
1398 istack[jstack - 1] = i;
1403 istack[jstack] = j - 1;
1404 istack[jstack - 1] = l;