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->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
98 pGrid =
new CSG_Grid(pGrid, Type, bCached);
100 if( !pGrid->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
108 if( !pGrid->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
114 CSG_Grid *pGrid =
new CSG_Grid(Type, NX, NY, Cellsize, xMin, yMin, bCached);
116 if( !pGrid->
is_Valid() ) {
delete(pGrid);
return( NULL ); }
return( pGrid );
203 void CSG_Grid::_On_Construction(
void)
209 m_Cache_Stream = NULL;
211 m_Cache_bSwap =
false;
212 m_Cache_bFlip =
false;
241 #pragma omp parallel for
242 for(
int y=0; y<
Get_NY(); y++)
244 for(
int x=0; x<
Get_NX(); x++)
296 if( _Load_PGSQL (File, bCached, bLoadData)
297 || _Load_Native (File, bCached, bLoadData)
298 || _Load_Compressed(File, bCached, bLoadData)
299 || _Load_Surfer (File, bCached, bLoadData)
300 || _Load_External (File, bCached, bLoadData) )
326 _Set_Properties(Type, NX, NY, Cellsize, xMin, yMin);
328 return( _Memory_Create(bCached) );
361 m_System .
Assign(0., 0., 0., 0, 0);
374 void CSG_Grid::_Set_Properties(
TSG_Data_Type Type,
int NX,
int NY,
double Cellsize,
double xMin,
double yMin)
397 m_System.
Assign(Cellsize > 0. ? Cellsize : 1., xMin, yMin, NX, NY);
416 if( (Scale != m_zScale && Scale != 0.) || Offset != m_zOffset )
453 return( m_Values != NULL ||
is_Cached() );
483 return( m_System == System );
503 if( y >= 0 && y <
Get_NY() )
507 for(
int x=0; x<
Get_NX(); x++)
521 for(
int x=0; x<
Get_NX(); x++)
556 return(
Get_Value(p.
x, p.
y, Value, Resampling, bNoData, bByteWise) );
567 if( bNoData ||
is_InGrid(ix + (
int)(0.5 + dx), iy + (
int)(0.5 + dy)) )
590 inline bool CSG_Grid::_Get_ValAtPos_NearestNeighbour(
double &Value,
int x,
int y,
double dx,
double dy)
const
592 if(
is_InGrid(x = x + (
int)(0.5 + dx), y = y + (
int)(0.5 + dy)) )
608 #define BILINEAR_ADD(ix, iy, d) if( is_InGrid(ix, iy) ) {\
609 n += d; z += d * asDouble(ix, iy);\
613 #define BILINEAR_ADD_BYTE(ix, iy, d) if( is_InGrid(ix, iy) ) {\
614 n += d; sLong v = asInt(ix, iy);\
615 z[0] += d * SG_GET_BYTE_0(v);\
616 z[1] += d * SG_GET_BYTE_1(v);\
617 z[2] += d * SG_GET_BYTE_2(v);\
618 z[3] += d * SG_GET_BYTE_3(v);\
622 inline bool CSG_Grid::_Get_ValAtPos_BiLinear(
double &Value,
int x,
int y,
double dx,
double dy,
bool bByteWise)
const
626 double z = 0., n = 0.;
671 inline double CSG_Grid::_Get_ValAtPos_BiCubicSpline(
double dx,
double dy,
double v_xy[4][4])
const
673 #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]))))
677 for(
int ix=0; ix<4; ix++)
686 inline bool CSG_Grid::_Get_ValAtPos_BiCubicSpline(
double &Value,
int x,
int y,
double dx,
double dy,
bool bByteWise)
const
692 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
694 Value = _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy);
701 double v_xy[4][4][4];
703 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
706 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[0]),
707 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[1]),
708 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[2]),
709 _Get_ValAtPos_BiCubicSpline(dx, dy, v_xy[3])
725 inline double CSG_Grid::_Get_ValAtPos_BSpline(
double dx,
double dy,
double v_xy[4][4])
const
729 for(
int i=0; i<4; i++)
734 if( (d = i - dx + 1.) > 0. ) s += d*d*d;
735 if( (d = i - dx + 0.) > 0. ) s += -4. * d*d*d;
736 if( (d = i - dx - 1.) > 0. ) s += 6. * d*d*d;
737 if( (d = i - dx - 2.) > 0. ) s += -4. * d*d*d;
742 if( (d = i - dy + 1.) > 0. ) s += d*d*d;
743 if( (d = i - dy + 0.) > 0. ) s += -4. * d*d*d;
744 if( (d = i - dy - 1.) > 0. ) s += 6. * d*d*d;
745 if( (d = i - dy - 2.) > 0. ) s += -4. * d*d*d;
752 for(
int iy=0; iy<4; iy++)
754 for(
int ix=0; ix<4; ix++)
756 z += v_xy[ix][iy] * Rx[ix] * Ry[iy];
764 inline bool CSG_Grid::_Get_ValAtPos_BSpline(
double &Value,
int x,
int y,
double dx,
double dy,
bool bByteWise)
const
770 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
772 Value = _Get_ValAtPos_BSpline(dx, dy, v_xy);
779 double v_xy[4][4][4];
781 if( _Get_ValAtPos_Fill4x4Submatrix(x, y, v_xy) )
784 _Get_ValAtPos_BSpline(dx, dy, v_xy[0]),
785 _Get_ValAtPos_BSpline(dx, dy, v_xy[1]),
786 _Get_ValAtPos_BSpline(dx, dy, v_xy[2]),
787 _Get_ValAtPos_BSpline(dx, dy, v_xy[3])
803 inline bool CSG_Grid::_Get_ValAtPos_Fill4x4Submatrix(
int x,
int y,
double v_xy[4][4])
const
805 int ix, iy, jx, jy, nNoData = 0;
808 for(iy=0, jy=y-1; iy<4; iy++, jy++)
810 for(ix=0, jx=x-1; ix<4; ix++, jx++)
826 for(
int i=0; nNoData>0 && nNoData<16 && i<16; i++)
830 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
832 t_xy[ix][iy] = v_xy[ix][iy];
835 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
842 for(jy=iy-1; jy<=iy+1; jy++)
for(jx=ix-1; jx<=ix+1; jx++)
846 s +=
asDouble(jx + x - 1, jy + y - 1);
849 else if( jy >= 0 && jy < 4 && jx >= 0 && jx < 4 && !
is_NoData_Value(t_xy[jx][jy]) )
858 v_xy[ix][iy] = s / n;
867 return( nNoData == 0 );
871 inline bool CSG_Grid::_Get_ValAtPos_Fill4x4Submatrix(
int x,
int y,
double v_xy[4][4][4])
const
873 int ix, iy, jx, jy, nNoData = 0;
876 for(iy=0, jy=y-1; iy<4; iy++, jy++)
878 for(ix=0, jx=x-1; ix<4; ix++, jx++)
882 int z =
asInt(jx, jy);
891 v_xy[0][ix][iy] = -1;
899 for(
int i=0; nNoData>0 && nNoData<16 && i<16; i++)
901 double t_xy[4][4][4];
903 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
905 t_xy[0][ix][iy] = v_xy[0][ix][iy];
906 t_xy[1][ix][iy] = v_xy[1][ix][iy];
907 t_xy[2][ix][iy] = v_xy[2][ix][iy];
908 t_xy[3][ix][iy] = v_xy[3][ix][iy];
911 for(iy=0; iy<4; iy++)
for(ix=0; ix<4; ix++)
913 if( t_xy[0][ix][iy] < 0 )
916 double s[4]; s[0] = s[1] = s[2] = s[3] = 0.;
918 for(jy=iy-1; jy<=iy+1; jy++)
for(jx=ix-1; jx<=ix+1; jx++)
922 int z =
asInt(jx + x - 1, jy + y - 1);
930 else if( jy >= 0 && jy < 4 && jx >= 0 && jx < 4 && !
is_NoData_Value(t_xy[0][jx][jy]) )
932 s[0] += t_xy[0][jx][jy];
933 s[1] += t_xy[1][jx][jy];
934 s[2] += t_xy[2][jx][jy];
935 s[3] += t_xy[3][jx][jy];
942 v_xy[0][ix][iy] = s[0] / n;
943 v_xy[1][ix][iy] = s[1] / n;
944 v_xy[2][ix][iy] = s[2] / n;
945 v_xy[3][ix][iy] = s[3] / n;
954 return( nNoData == 0 );
989 m_Statistics += Scaling ? Offset + Scaling * Value : Value;
1001 for(
int x=0; x<
Get_NX(); x++)
1003 double Value =
asDouble(x, y,
false);
1007 m_Statistics += Scaling ? Offset + Scaling * Value : Value;
1063 if( Quantile <= 0. ) {
return(
Get_Min() ); }
1064 if( Quantile >= 1. ) {
return(
Get_Max() ); }
1066 if( bFromHistogram )
1086 return(
Get_Quantile(0.01 * Percentile, bFromHistogram) );
1103 Update();
return( m_Statistics );
1120 if( xMin > xMax || yMin > yMax )
1125 Statistics.
Create(bHoldValues);
1127 int nx = 1 + (xMax - xMin);
1128 int ny = 1 + (yMax - yMin);
1129 sLong nCells = nx * ny;
1137 for(
double i=0; i<(double)nCells; i+=d)
1139 int y = yMin + (int)i / nx;
1140 int x = xMin + (int)i % nx;
1142 double Value =
asDouble(x, y,
false);
1146 Statistics += Scaling ? Offset + Scaling * Value : Value;
1152 for(
int y=yMin; y<=yMax; y++)
1154 for(
int x=xMin; x<=xMax; x++)
1156 double Value =
asDouble(x, y,
false);
1160 Statistics += Scaling ? Offset + Scaling * Value : Value;
1170 #define SG_GRID_HISTOGRAM_CLASSES_DEFAULT 255
1191 return( m_Histogram );
1209 if( xMin > xMax || yMin > yMax )
1216 int nx = 1 + (xMax - xMin);
1217 int ny = 1 + (yMax - yMin);
1218 sLong nCells = nx * ny;
1226 for(
double i=0; i<(double)nCells; i+=d)
1228 int y = yMin + (int)i / nx;
1229 int x = xMin + (int)i % nx;
1231 double Value =
asDouble(x, y,
false);
1235 Histogram += Scaling ? Offset + Scaling * Value : Value;
1241 for(
int y=yMin; y<=yMax; y++)
1243 for(
int x=xMin; x<=xMax; x++)
1245 double Value =
asDouble(x, y,
false);
1249 Histogram += Scaling ? Offset + Scaling * Value : Value;
1255 return( Histogram.
Update() );
1266 bool CSG_Grid::_Set_Index(
void)
1279 sLong i, j, k, l, ir, n, *istack, jstack, nstack, indxt, itemp, nData;
1289 m_Index[--nData] = i;
1328 for(j=l+1; j<=ir; j++)
1333 for(i=j-1; i>=0; i--)
1340 m_Index[i + 1] = m_Index[i];
1343 m_Index[i + 1] = indxt;
1351 ir = istack[jstack--];
1352 l = istack[jstack--];
1358 #define SORT_SWAP(a,b) {itemp=(a);(a)=(b);(b)=itemp;}
1380 do i++;
while(
asDouble(m_Index[i]) < a);
1381 do j--;
while(
asDouble(m_Index[j]) > a);
1393 m_Index[l] = m_Index[j];
1397 if( jstack >= nstack )
1403 if( ir - i + 1 >= j - l )
1405 istack[jstack] = ir;
1406 istack[jstack - 1] = i;
1411 istack[jstack] = j - 1;
1412 istack[jstack - 1] = l;