Statistics
| Revision:

gvsig-raster / org.gvsig.raster / trunk / org.gvsig.raster / org.gvsig.raster.lib / org.gvsig.raster.lib.impl / src / main / java / org / gvsig / raster / impl / grid / GridInterpolated.java @ 2438

History | View | Annotate | Download (14.2 KB)

1
/*******************************************************************************
2
                GridWrapperInterpolated
3
                Copyright (C) Victor Olaya
4
                
5
                This program is free software; you can redistribute it and/or modify
6
                it under the terms of the GNU General Public License as published by
7
                the Free Software Foundation; either version 2 of the License, or
8
                (at your option) any later version.
9

10
                This program is distributed in the hope that it will be useful,
11
                but WITHOUT ANY WARRANTY; without even the implied warranty of
12
                MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
                GNU General Public License for more details.
14

15
                You should have received a copy of the GNU General Public License
16
                along with this program; if not, write to the Free Software
17
                Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
18
*******************************************************************************/ 
19
package org.gvsig.raster.impl.grid;
20

    
21
import java.awt.Rectangle;
22

    
23
import org.gvsig.fmap.dal.coverage.dataset.Buffer;
24
import org.gvsig.fmap.dal.coverage.datastruct.GridExtent;
25
import org.gvsig.fmap.dal.coverage.exception.ProcessInterruptedException;
26
import org.gvsig.fmap.dal.coverage.exception.QueryException;
27
import org.gvsig.fmap.dal.coverage.exception.RasterBufferInvalidAccessException;
28
import org.gvsig.fmap.dal.coverage.exception.RasterBufferInvalidException;
29
import org.gvsig.fmap.dal.coverage.store.RasterDataStore;
30
import org.gvsig.fmap.dal.coverage.store.RasterQuery;
31
import org.gvsig.raster.impl.DefaultRasterManager;
32

    
33
/**
34
 * A grid wrapper that performs interpolation to calculate
35
 * cell values. This should be used when the window extent
36
 * does not 'fit' into the structure (coordinates and cellsize)
37
 * of the grid.
38
 * 
39
 * @author Victor Olaya
40
 */
41

    
42
public class GridInterpolated extends GridReader {
43

    
44
        public static final int INTERPOLATION_NearestNeighbour = Buffer.INTERPOLATION_NearestNeighbour;
45
        public static final int INTERPOLATION_Bilinear = Buffer.INTERPOLATION_Bilinear;
46
        public static final int INTERPOLATION_InverseDistance = Buffer.INTERPOLATION_InverseDistance;
47
        public static final int INTERPOLATION_BicubicSpline = Buffer.INTERPOLATION_BicubicSpline;
48
        public static final int INTERPOLATION_BSpline = Buffer.INTERPOLATION_BSpline;
49
        
50
        double m_dXMin , m_dYMax; //coordinates of the layer, not the window
51
        double m_dCellSizeX; //cellsize of the layer, not the window
52
        double m_dCellSizeY; //cellsize of the layer, not the window
53
        int m_iInterpolationMethod = INTERPOLATION_BSpline;
54

    
55
        /**
56
         * Crea un objeto lector a partir de un buffer de datos y el extent de la extensi?n
57
         * completa y de la ventana accedida.
58
         * @param rb Buffer de datos
59
         * @param layerExtent extent de la capa completa
60
         * @param windowExtent Extent
61
         * @param bands N?mero de bandas del origen
62
         */
63
        public GridInterpolated(        Buffer rb, 
64
                                                                GridExtent layerExtent,
65
                                                                GridExtent windowExtent,
66
                                                                int[] bands){
67
                super(rb, layerExtent, windowExtent, bands);
68
                m_dXMin = layerExtent.minX();
69
                m_dYMax = layerExtent.maxY();
70
                
71
                m_dCellSizeX = layerExtent.getCellSizeX();
72
                m_dCellSizeY = layerExtent.getCellSizeY();
73
        }
74
        
75
        /**
76
         * Crea un objeto lector a partir de una fuente de datos y el extent de la extensi?n
77
         * completa y de la ventana accedida.
78
         * @param ds Fuente de datos
79
         * @param layerExtent extent de la capa completa
80
         * @param windowExtent Extent
81
         * @param bands N?mero de bandas del origen
82
         */
83
        public GridInterpolated(        RasterDataStore ds, 
84
                                                                GridExtent layerExtent,
85
                                                                GridExtent windowExtent,
86
                                                                int[] bands){
87
                super(ds, layerExtent, windowExtent, bands);
88
                init();
89
        }
90
        
91
        /**
92
         * Inicializaci?n de la fuente de datos y carga del buffer.
93
         */
94
        private void init(){
95
                m_dXMin = dataStore.getExtent().minX();
96
                m_dYMax = dataStore.getExtent().maxY();
97
                
98
                m_dCellSizeX = dataStore.getAffineTransform().getScaleX();
99
                m_dCellSizeY = dataStore.getAffineTransform().getScaleY();
100
                
101
        }
102
                
103
        /**
104
         * Asigna el m?todo de interpolaci?n 
105
         * @param iMethod
106
         */
107
        public void setInterpolationMethod(int iMethod){
108
                m_iInterpolationMethod = iMethod;
109
        }
110

    
111
        /*
112
         *  (non-Javadoc)
113
         * @see org.gvsig.fmap.grid.GridReader#getCellValueAsByte(int, int)
114
         */
115
        public byte getCellValueAsByte(int x, int y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
116
                return  (byte)_getValueAt(x, y);
117
        }
118
        
119
        /*
120
         *  (non-Javadoc)
121
         * @see org.gvsig.fmap.grid.GridReader#getCellValueAsShort(int, int)
122
         */
123
        public short getCellValueAsShort(int x, int y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
124
                return (short) _getValueAt(x, y);
125
        }
126
        
127
        /*
128
         *  (non-Javadoc)
129
         * @see org.gvsig.fmap.grid.GridReader#getCellValueAsInt(int, int)
130
         */
131
        public int getCellValueAsInt(int x, int y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {                
132
                return (int) _getValueAt(x, y);
133
        }
134

    
135
        /*
136
         *  (non-Javadoc)
137
         * @see org.gvsig.fmap.grid.GridReader#getCellValueAsFloat(int, int)
138
         */
139
        public float getCellValueAsFloat(int x, int y)  throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
140
                return (float) _getValueAt(x, y);
141
        }
142
        
143
        /*
144
         *  (non-Javadoc)
145
         * @see org.gvsig.fmap.grid.GridReader#getCellValueAsDouble(int, int)
146
         */
147
        public double getCellValueAsDouble(int x, int y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
148
                return _getValueAt(x, y);
149
        }
150
        
151
        /*
152
         *  (non-Javadoc)
153
         * @see org.gvsig.fmap.grid.GridReader#getCellValueAsDouble(int, int)
154
         */
155
        public double getCellValue(int x, int y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
156
                return _getValueAt(x, y);
157
        }
158
        
159
        private double _getValueAt(int x, int y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
160
                double dX = windowExtent.minX() + windowExtent.getCellSize() * (x + 0.5);
161
                double dY = windowExtent.maxY() - windowExtent.getCellSize() * (y + 0.5);
162
                return getValueAt(dX, dY);
163
        }
164

    
165
        
166
        /** Devuelve el valor interpolado para unas coordenadas de la imagen no enteras
167
         *  la distancia se toma al vertice superior izquierdo del pixel
168
         * */
169
        public double _getValueAt(double x, double y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
170
                int x_int, y_int;
171
                x_int= (int)Math.floor(x);
172
                y_int= (int)Math.floor(y);
173
                double dX = windowExtent.minX() + windowExtent.getCellSize() * (x_int + (x-x_int));
174
                double dY = windowExtent.maxY() - windowExtent.getCellSize() * (y_int + (y-y_int));
175
                return getValueAt(dX, dY);
176
        }
177
        
178
        
179
        private double getValueAt(double xPosition, double yPosition) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
180
                int        x, y;
181
                double        dx, dy;
182
                double dValue;
183

    
184
                x        = (int) Math.floor(xPosition = (xPosition - m_dXMin) / m_dCellSizeX);
185
                y        = (int) Math.floor(yPosition = (yPosition - m_dYMax) / m_dCellSizeY);
186

    
187
                dValue = _getCellValueInLayerCoords(x, y);
188
                
189
                if(        !isNoDataValue(dValue) ) {
190

    
191
                        dx        = xPosition - x;
192
                        dy        = yPosition - y;
193
                        
194
                        switch( m_iInterpolationMethod ) {
195
                        case INTERPOLATION_NearestNeighbour:
196
                                dValue = _getValueNearestNeighbour (x, y, dx, dy);
197
                                break;
198

    
199
                        case INTERPOLATION_Bilinear:
200
                                dValue        = _getValueBiLinear (x, y, dx, dy);
201
                                break;
202

    
203
                        case INTERPOLATION_InverseDistance:
204
                                dValue        = _getValueInverseDistance(x, y, dx, dy);
205
                                break;
206

    
207
                        case INTERPOLATION_BicubicSpline:
208
                                dValue        = _getValueBiCubicSpline (x, y, dx, dy);
209
                                break;
210

    
211
                        case INTERPOLATION_BSpline:
212
                                dValue        = _getValueBSpline (x, y, dx, dy);
213
                                break;
214
                        }
215
                } else {
216
                        dValue = getNoDataValue();
217
                }
218

    
219
                return dValue;
220
        }
221

    
222
        private double _getValueNearestNeighbour(int x, int y, double dx, double dy) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
223
                x += (int)(0.5 + dx);
224
                y += (int)(0.5 + dy);
225

    
226
                return _getCellValueInLayerCoords(x, y);
227
        }
228

    
229
        private  double _getValueBiLinear(int x, int y, double dx, double dy) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
230
                double        z = 0.0, n = 0.0, d;
231
                double dValue;
232
                
233
                dValue = _getCellValueInLayerCoords(x, y);
234
                if (!isNoDataValue(dValue)){
235
                         d = (1.0 - dx) * (1.0 - dy);
236
                         z += d * dValue; 
237
                         n += d;
238
                }
239
                
240
                dValue = _getCellValueInLayerCoords(x + 1, y);
241
                if (!isNoDataValue(dValue)){
242
                         d = (dx) * (1.0 - dy);
243
                         z += d * dValue; 
244
                         n += d;
245
                }
246
                
247
                dValue = _getCellValueInLayerCoords(x, y + 1);
248
                if (!isNoDataValue(dValue)){
249
                         d = (1.0 - dx) * (dy);
250
                         z += d * dValue; 
251
                         n += d;
252
                }
253
                
254
                dValue = _getCellValueInLayerCoords(x + 1, y + 1);
255
                if (!isNoDataValue(dValue)){
256
                         d = (dx) * (dy);
257
                         z += d * dValue; 
258
                         n += d;
259
                }
260

    
261
                if( n > 0.0 ){
262
                        return( z / n );
263
                }
264

    
265
                return( getNoDataValue() );
266
        }
267

    
268
        private double _getValueInverseDistance(int x, int y, double dx, double dy) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
269
                double        z = 0.0, n = 0.0, d;
270
                double dValue;
271
                
272
                if( dx > 0.0 || dy > 0.0 ){
273

    
274
                        dValue = _getCellValueInLayerCoords(x, y);
275
                        if (!isNoDataValue(dValue)){
276
                                d = 1.0 / Math.sqrt(dx*dx + dy*dy); 
277
                                z += d * dValue; 
278
                                n += d;
279
                        }
280
                        
281
                        dValue = _getCellValueInLayerCoords(x + 1, y);
282
                        if (!isNoDataValue(dValue)){
283
                                d = 1.0 / Math.sqrt((1.0-dx)*(1.0-dx) + dy*dy); 
284
                                z += d * dValue; 
285
                                n += d;
286
                        }
287
                        
288
                        dValue = _getCellValueInLayerCoords(x, y + 1);
289
                        if (!isNoDataValue(dValue)){
290
                                d = 1.0 / Math.sqrt(dx*dx + (1.0-dy)*(1.0-dy)); 
291
                                z += d * dValue; 
292
                                n += d;
293
                        }
294
                        
295
                        dValue = _getCellValueInLayerCoords(x + 1, y + 1);
296
                        if (!isNoDataValue(dValue)){
297
                                d = 1.0 / Math.sqrt((1.0-dx)*(1.0-dx) + (1.0-dy)*(1.0-dy)); 
298
                                z += d * dValue; 
299
                                n += d;
300
                        }
301

    
302
                        if( n > 0.0 ){
303
                                return( z / n );
304
                        }
305
                }else
306
                        return( _getCellValueInLayerCoords(x, y) );
307
                
308
                return( getNoDataValue());
309
        }
310

    
311
        private double _getValueBiCubicSpline(int x, int y, double dx, double dy) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
312
                int                i;
313
                double        a0, a2, a3, b1, b2, b3, c[], z_xy[][];
314

    
315
                c = new double[4];
316
                z_xy = new double[4][4];
317
                
318
                if( _get4x4Submatrix(x, y, z_xy) ){
319
                        
320
                        for(i=0; i<4; i++){
321
                                a0                = z_xy[0][i] - z_xy[1][i];
322
                                a2                = z_xy[2][i] - z_xy[1][i];
323
                                a3                = z_xy[3][i] - z_xy[1][i];
324

    
325
                                b1                = -a0 / 3.0 + a2       - a3 / 6.0;
326
                                b2                =  a0 / 2.0 + a2 / 2.0;
327
                                b3                = -a0 / 6.0 - a2 / 2.0 + a3 / 6.0;
328

    
329
                                c[i]        = z_xy[1][i] + b1 * dx + b2 * dx*dx + b3 * dx*dx*dx;
330
                        }
331

    
332
                        a0                = c[0] - c[1];
333
                        a2                = c[2] - c[1];
334
                        a3                = c[3] - c[1];
335

    
336
                        b1                = -a0 / 3.0 + a2       - a3 / 6.0;
337
                        b2                =  a0 / 2.0 + a2 / 2.0;
338
                        b3                = -a0 / 6.0 - a2 / 2.0 + a3 / 6.0;
339

    
340
                        return( c[1] + b1 * dy + b2 * dy*dy + b3 * dy*dy*dy );
341
                }
342

    
343
                return( _getValueBiLinear(x, y, dx, dy) );
344
        }
345

    
346
        private double _getValueBSpline(int x, int y, double dx, double dy) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
347
                int                i, ix, iy;
348
                double        z, px, py, Rx[], Ry[], z_xy[][];
349
                
350
                Rx = new double[4];
351
                Ry = new double[4];
352
                z_xy = new double [4][4];
353

    
354
                if( _get4x4Submatrix(x, y, z_xy) ) {
355
                        for(i=0, px=-1.0-dx, py=-1.0-dy; i<4; i++, px++, py++) {
356
                                Rx[i]        = 0.0;
357
                                Ry[i]        = 0.0;
358

    
359
                                if( (z = px + 2.0) > 0.0 )
360
                                        Rx[i]        +=        z*z*z;
361
                                if( (z = px + 1.0) > 0.0 )
362
                                        Rx[i]        += -4.0 * z*z*z;
363
                                if( (z = px + 0.0) > 0.0 )
364
                                        Rx[i]        +=  6.0 * z*z*z;
365
                                if( (z = px - 1.0) > 0.0 )
366
                                        Rx[i]        += -4.0 * z*z*z;
367
                                if( (z = py + 2.0) > 0.0 )
368
                                        Ry[i]        +=        z*z*z;
369
                                if( (z = py + 1.0) > 0.0 )
370
                                        Ry[i]        += -4.0 * z*z*z;
371
                                if( (z = py + 0.0) > 0.0 )
372
                                        Ry[i]        +=  6.0 * z*z*z;
373
                                if( (z = py - 1.0) > 0.0 )
374
                                        Ry[i]        += -4.0 * z*z*z;
375

    
376
                                Rx[i]        /= 6.0;
377
                                Ry[i]        /= 6.0;
378
                        }
379

    
380
                        for(iy=0, z=0.0; iy<4; iy++) {
381
                                for(ix=0; ix<4; ix++){
382
                                        z        += z_xy[ix][iy] * Rx[ix] * Ry[iy];
383
                                }
384
                        }
385
                        return( z );
386
                }
387
                return( _getValueBiLinear(x, y, dx, dy) );
388
        }
389

    
390
        private boolean _get4x4Submatrix(int x, int y, double z_xy[][]) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
391
                int        ix, iy, px, py;
392
                double dValue;
393

    
394
                for(iy=0, py=y-1; iy<4; iy++, py++) {
395
                        for(ix=0, px=x-1; ix<4; ix++, px++) {
396
                                dValue = _getCellValueInLayerCoords(px, py);
397
                                if (isNoDataValue(dValue)) {
398
                                        return false;
399
                                } else {
400
                                        z_xy[ix][iy] = dValue;
401
                                }
402
                        }
403
                }
404
                return( true );
405
        }
406

    
407
        private double _getCellValueInLayerCoords(int x, int y) throws RasterBufferInvalidAccessException, RasterBufferInvalidException {
408
                int w, h;
409
                int pX = 0, pY = 0;
410
                RasterQuery query = DefaultRasterManager.getInstance().createQuery();
411
                query.setDrawableBands(bands);
412
                query.storeLastBuffer(true);
413
                
414
                if (dataStore != null) {
415
                        w = (int)dataStore.getWidth();
416
                        h = (int)dataStore.getHeight();
417
                        try {
418
                                if (y >= h || x >= w || x < 0 || y < 0)
419
                                        return getNoDataValue();
420
                                query.setAreaOfInterest(new Rectangle(x, y, 1, 1));
421
                                rasterBuf = dataStore.query(query);
422
                        } catch (QueryException e) {
423
                                throw new RasterBufferInvalidAccessException("");
424
                        } catch (ProcessInterruptedException e) {
425
                                // La cancelaci?n de la lectura de un pixel no es significativa
426
                        }
427
                } else {
428
                        w = rasterBuf.getWidth();
429
                        h = rasterBuf.getHeight();
430
                        pX = x;
431
                        pY = y;
432
                }
433

    
434
                if (x >= 0 && y >= 0 && x < w && y < h) {
435
                        switch (dataType) {
436
                                case Buffer.TYPE_DOUBLE:
437
                                        return rasterBuf.getElemDouble(pY, pX, bandToOperate);
438
                                case Buffer.TYPE_INT:
439
                                        return (double) rasterBuf.getElemInt(pY, pX, bandToOperate);
440
                                case Buffer.TYPE_FLOAT:
441
                                        return (double) rasterBuf.getElemFloat(pY, pX, bandToOperate);
442
                                case Buffer.TYPE_BYTE:
443
                                        return ((byte) (rasterBuf.getElemByte(pY, pX, bandToOperate)) & 0xff);
444
                                case Buffer.TYPE_SHORT:
445
                                case Buffer.TYPE_USHORT:
446
                                        return (double) rasterBuf.getElemShort(pY, pX, bandToOperate);
447
                        }
448
                }
449
                return getNoDataValue();
450
        }
451

    
452
        public byte[] getBandsValuesAsByte(int x, int y) {
453
                //TODO: FUNCIONALIDAD: getBandsValuesAsByte:Obtener los valores interpolados en todas las bandas
454
                return null;
455
        }
456

    
457
        public short[] getBandsValuesAsShort(int x, int y) {
458
                //TODO: FUNCIONALIDAD: getBandsValuesAsShort:Obtener los valores interpolados en todas las bandas
459
                return null;
460
        }
461

    
462
        public int[] getBandsValuesAsInt(int x, int y) {
463
                //TODO: FUNCIONALIDAD: getBandsValuesAsInt:Obtener los valores interpolados en todas las bandas
464
                return null;
465
        }
466

    
467
        public float[] getBandsValuesAsFloat(int x, int y) {
468
                //TODO: FUNCIONALIDAD: getBandsValuesAsFloat:Obtener los valores interpolados en todas las bandas
469
                return null;
470
        }
471

    
472
        public double[] getBandsValuesAsDouble(int x, int y) {
473
                //TODO: FUNCIONALIDAD: getBandsValuesAsDouble:Obtener los valores interpolados en todas las bandas
474
                return null;
475
        }
476
}