gistools/readers/geotiff/
image.rs

1use super::decoder::Decoder;
2use crate::{
3    parsers::{RGBA, Reader},
4    proj::Transformer,
5    readers::{
6        FieldTagNames, GTiffDataType, GeoKeyDirectoryKeys, GeoTIFFVectorFeature, GeoTiePoint,
7        ImageDirectory, PhotometricInterpretations, apply_predictor, build_samples,
8        build_transform_from_geo_keys, convert_color_space, geotiff::decoder::get_decoder,
9        get_reader_for_sample, needs_normalization, normalize_array, sample_sum,
10    },
11};
12use alloc::{rc::Rc, vec, vec::Vec};
13use core::cell::RefCell;
14use half::f16;
15use libm::{fmax, fmin, pow, round, sqrt};
16use s2json::{BBox, Properties, VectorGeometry, VectorMultiPoint, VectorPoint};
17
18/// Metadata for a GeoTIFF image
19#[derive(Debug, Clone, Default, PartialEq)]
20pub struct GeoTIFFMetadata {
21    /// The height of the image
22    pub height: usize,
23    /// The width of the image
24    pub width: usize,
25    /// True if the image has an alpha channel
26    pub alpha: bool,
27}
28
29/// Raster data
30#[derive(Debug, Clone, Default, PartialEq)]
31pub struct Raster {
32    /// The width of the image
33    pub width: usize,
34    /// The height of the image
35    pub height: usize,
36    /// Type
37    pub r#type: GTiffDataType,
38    /// The type data
39    pub data: Vec<f64>,
40    /// True if the image has an alpha channel
41    pub alpha: bool,
42    /// The min value found
43    pub min: f64,
44    /// The max value found
45    pub max: f64,
46}
47impl Raster {
48    /// Convert the data to u8s
49    pub fn to_u8s(&self) -> Vec<u8> {
50        self.data.iter().map(|x| round(fmax(0., fmin(255., *x))) as u8).collect()
51    }
52    /// Convert the data to u16s
53    pub fn to_u16s(&self) -> Vec<u16> {
54        self.data.iter().map(|x| round(fmax(0., fmin(65535., *x))) as u16).collect()
55    }
56    /// Convert the data to u32s
57    pub fn to_u32s(&self) -> Vec<u32> {
58        self.data.iter().map(|x| round(fmax(0., fmin(4294967295., *x))) as u32).collect()
59    }
60    /// Convert the data to i8s
61    pub fn to_i8s(&self) -> Vec<i8> {
62        self.data.iter().map(|x| round(fmax(-128., fmin(127., *x))) as i8).collect()
63    }
64    /// Convert the data to i16s
65    pub fn to_i16s(&self) -> Vec<i16> {
66        self.data.iter().map(|x| round(fmax(-32768., fmin(32767., *x))) as i16).collect()
67    }
68    /// Convert the data to i32s
69    pub fn to_i32s(&self) -> Vec<i32> {
70        self.data.iter().map(|x| round(fmax(-2147483648., fmin(2147483647., *x))) as i32).collect()
71    }
72    /// Convert the data to f16s
73    pub fn to_f16s(&self) -> Vec<f16> {
74        self.data.iter().map(|x| f16::from_f64(*x)).collect()
75    }
76    /// Convert the data to f32s
77    pub fn to_f32s(&self) -> Vec<f32> {
78        self.data.iter().map(|x| *x as f32).collect()
79    }
80}
81
82/// Internal interface for sample reader
83pub type SampleReader = fn(buffer: &[u8], offset: usize, little_endian: bool) -> f64;
84
85/// Internal interface for sample format
86pub struct SampleFormat {
87    /// The source sample offset collection
88    pub src_sample_offsets: Vec<u32>,
89    /// The sample readers
90    pub sample_readers: Vec<SampleReader>,
91}
92
93/// # GeoTIFF Image Container
94///
95/// ## Description
96/// A Container for a GeoTIFF image
97///
98/// ## Usage
99///
100/// The methods you have access to:
101/// - [`GeoTIFFImage::new`]: Create a new GeoTIFFImage
102/// - [`GeoTIFFImage::width`]: Get the image width
103/// - [`GeoTIFFImage::height`]: Get the image height
104/// - [`GeoTIFFImage::tile_width`]: Get the tile width
105/// - [`GeoTIFFImage::tile_height`]: Get the tile height
106/// - [`GeoTIFFImage::block_width`]: Get the block width
107/// - [`GeoTIFFImage::block_height`]: Get the block height
108/// - [`GeoTIFFImage::bytes_per_pixel`]: Calculates the number of bytes for each pixel across all samples.
109/// - [`GeoTIFFImage::samples_per_pixel`]: The number of samples per pixel
110/// - [`GeoTIFFImage::get_sample_format`]: Returns the sample format
111/// - [`GeoTIFFImage::get_bits_per_sample`]: Returns the number of bits per sample
112/// - [`GeoTIFFImage::raster_array_type`]: Convert the data format and bits per sample to the appropriate array type
113/// - [`GeoTIFFImage::get_tie_point`]: Returns an array of tiepoints.
114/// - [`GeoTIFFImage::origin`]: Returns the image origin as a XYZ-vector.
115/// - [`GeoTIFFImage::origin_ll`]: Returns the image origin as a XYZ-vector in lon-lat space.
116/// - [`GeoTIFFImage::resolution`]: Returns the image resolution as a XYZ-vector.
117/// - [`GeoTIFFImage::resolution_ll`]: Returns the image resolution as a XYZ-vector in lon-lat space.
118/// - [`GeoTIFFImage::pixel_is_area`]: Returns whether or not the pixels of the image depict an area (or point).
119/// - [`GeoTIFFImage::get_bbox`]: Returns the image bounding box as an array of 4 values: min-x, min-y, max-x and max-y.
120/// - [`GeoTIFFImage::raster_data`]: Returns the raster data of the image.
121/// - [`GeoTIFFImage::get_rgba`]: Returns the RGBA raster data of the image.
122/// - [`GeoTIFFImage::get_multi_point_vector`]: Build a vector feature from the image
123///
124/// The methods you probably care about are `get_bbox`, `raster_data`, `get_rgba`,
125/// and `get_multi_point_vector`
126#[derive(Debug, Clone, Default)]
127pub struct GeoTIFFImage<T: Reader> {
128    reader: Rc<RefCell<T>>,
129    image_directory: ImageDirectory,
130    little_endian: bool,
131    is_tiled: bool,
132    transformer: Rc<RefCell<Transformer>>,
133    decode_fn: Option<Decoder>,
134    src_sample_offsets: Vec<u32>,
135    sample_readers: Vec<SampleReader>,
136    planar_configuration: i16,
137}
138impl<T: Reader> GeoTIFFImage<T> {
139    /// Create a new GeoTIFFImage
140    pub fn new(
141        reader: Rc<RefCell<T>>,
142        image_directory: ImageDirectory,
143        transformer: Rc<RefCell<Transformer>>,
144        little_endian: bool,
145    ) -> Self {
146        build_transform_from_geo_keys(
147            &mut transformer.borrow_mut(),
148            &image_directory.geo_key_directory,
149        );
150        let compression = image_directory.variables.get_short(FieldTagNames::Compression as u16);
151        let planar_configuration = image_directory
152            .variables
153            .get_short(FieldTagNames::PlanarConfiguration as u16)
154            .unwrap_or(1);
155        let is_tiled = !image_directory.variables.has(FieldTagNames::StripOffsets as u16);
156        Self {
157            reader,
158            image_directory,
159            little_endian,
160            is_tiled,
161            transformer,
162            decode_fn: get_decoder(compression.map(|v| v as u16)),
163            src_sample_offsets: Vec::new(),
164            sample_readers: Vec::new(),
165            planar_configuration,
166        }
167    }
168
169    /// Get the image width
170    pub fn width(&self) -> usize {
171        self.image_directory.variables.get_short(FieldTagNames::ImageWidth as u16).unwrap_or(0)
172            as usize
173    }
174
175    /// Get the image height
176    pub fn height(&self) -> usize {
177        self.image_directory.variables.get_short(FieldTagNames::ImageLength as u16).unwrap_or(0)
178            as usize
179    }
180
181    /// Get the tile width
182    pub fn tile_width(&self) -> usize {
183        if self.is_tiled {
184            self.image_directory.variables.get_short(FieldTagNames::TileWidth as u16).unwrap_or(0)
185                as usize
186        } else {
187            self.width()
188        }
189    }
190
191    /// Get the tile height
192    pub fn tile_height(&self) -> usize {
193        if self.is_tiled {
194            self.image_directory.variables.get_short(FieldTagNames::TileLength as u16).unwrap_or(0)
195                as usize
196        } else {
197            let rows_per_strip = self
198                .image_directory
199                .variables
200                .get_short(FieldTagNames::RowsPerStrip as u16)
201                .unwrap_or(0) as usize;
202            self.height().min(rows_per_strip)
203        }
204    }
205
206    /// Get the block width
207    pub fn block_width(&self) -> usize {
208        self.tile_width()
209    }
210
211    /// Get the block height
212    pub fn block_height(&self, y: usize) -> usize {
213        let tile_height = self.tile_height();
214        let height = self.height();
215        if self.is_tiled || (y + 1) * tile_height <= height {
216            tile_height
217        } else {
218            height - y * tile_height
219        }
220    }
221
222    /// Calculates the number of bytes for each pixel across all samples. Only full
223    /// bytes are supported, an exception is thrown when this is not the case.
224    ///
225    /// ## Returns
226    /// the bytes per pixel
227    pub fn bytes_per_pixel(&self) -> usize {
228        let bits_per_sample = self
229            .image_directory
230            .variables
231            .get_u16s(FieldTagNames::BitsPerSample as u16)
232            .unwrap_or_default();
233        let mut bytes = 0;
234        for bit in &bits_per_sample {
235            bytes += (*bit).div_ceil(8) as usize;
236        }
237        bytes
238    }
239
240    /// ## Returns
241    /// The number of samples per pixel
242    pub fn samples_per_pixel(&self) -> usize {
243        self.image_directory.variables.get_short(FieldTagNames::SamplesPerPixel as u16).unwrap_or(1)
244            as usize
245    }
246
247    /// Returns the sample format
248    ///
249    /// ## Parameters
250    /// - `sample_index`: the sample index to start at
251    ///
252    /// ## Returns
253    /// The sample format code
254    pub fn get_sample_format(&self, sample_index: Option<usize>) -> u16 {
255        let sample_index = sample_index.unwrap_or(0);
256        *self
257            .image_directory
258            .variables
259            .get_u16s(FieldTagNames::SampleFormat as u16)
260            .unwrap_or_default()
261            .get(sample_index)
262            .unwrap_or(&1)
263    }
264
265    /// Returns the number of bits per sample
266    ///
267    /// ## Parameters
268    /// - `sample_index`: the sample index to start at
269    ///
270    /// ## Returns
271    /// The number of bits per sample at the sample index
272    pub fn get_bits_per_sample(&self, sample_index: Option<usize>) -> u16 {
273        let sample_index = sample_index.unwrap_or(0);
274        *self
275            .image_directory
276            .variables
277            .get_u16s(FieldTagNames::BitsPerSample as u16)
278            .unwrap_or_default()
279            .get(sample_index)
280            .unwrap_or(&0)
281    }
282
283    /// Convert the data format and bits per sample to the appropriate array type
284    ///
285    /// ## Parameters
286    /// - `raster`: the data
287    ///
288    /// ## Returns
289    /// The array
290    pub fn raster_array_type(&self) -> GTiffDataType {
291        let format = self.get_sample_format(None);
292        let bits_per_sample = self.get_bits_per_sample(None);
293        GTiffDataType::to_type(format, bits_per_sample)
294    }
295
296    /// Returns an array of tiepoints.
297    pub fn get_tie_point(&self) -> GeoTiePoint {
298        self.image_directory.tie_point
299    }
300
301    /// Returns the image origin as a XYZ-vector. When the image has no affine
302    /// transformation, then an exception is thrown.
303    ///
304    /// ## Returns
305    /// The origin as a vector
306    pub fn origin(&self) -> VectorPoint<()> {
307        // const { tiepoint, ModelTransformation: transform } = self.#imageDirectory;
308        // if (Array.isArray(tiepoint) && tiepoint.length === 6) {
309        //   return { x: tiepoint[3], y: tiepoint[4], z: tiepoint[5] };
310        // } else if (transform !== undefined) {
311        //   return { x: transform[3], y: transform[7], z: transform[11] };
312        // }
313        // throw new Error('The image does not have an affine transformation.');
314        // TODO: Is a modeltranasformation needed here?
315        let tiepoint = self.get_tie_point();
316        VectorPoint::new_xyz(tiepoint.x, tiepoint.y, tiepoint.z, None)
317    }
318
319    /// Returns the image origin as a XYZ-vector in lon-lat space. When the image has no affine
320    /// transformation, then an exception is thrown.
321    ///
322    /// ## Returns
323    /// The origin as a lon-lat vector
324    pub fn origin_ll(&self) -> VectorPoint<()> {
325        let mut origin = self.origin();
326        self.transformer.borrow_mut().forward_mut(&mut origin);
327        origin
328    }
329
330    /// Returns the image resolution as a XYZ-vector. When the image has no affine
331    /// transformation, then an exception is thrown. in cases when the current image does
332    /// not have the required tags on its own.
333    ///
334    /// ## Returns
335    /// The resolution as a vector
336    pub fn resolution(&self) -> VectorPoint<()> {
337        // try pixel scale first
338        let pixel_scale = self.image_directory.pixel_scale;
339        if pixel_scale.x != 0. || pixel_scale.y != 0. || pixel_scale.z != 0. {
340            return VectorPoint::new_xyz(pixel_scale.x, -pixel_scale.y, pixel_scale.z, None);
341        }
342        // then try model transformation
343        let transform = self
344            .image_directory
345            .variables
346            .getf64s(FieldTagNames::ModelTransformation as u16)
347            .unwrap_or_else(|| panic!("The image does not have an affine transformation."));
348        if transform[1] == 0. && transform[4] == 0. {
349            VectorPoint::new_xyz(transform[0], -transform[5], transform[10], None)
350        } else {
351            let x = sqrt(transform[0] * transform[0] + transform[4] * transform[4]);
352            let y = -sqrt(transform[1] * transform[1] + transform[5] * transform[5]);
353            let z = transform[10];
354            VectorPoint::new_xyz(x, y, z, None)
355        }
356    }
357
358    /// Returns the image resolution as a XYZ-vector in lon-lat space. When the image has no affine
359    /// transformation, then an exception is thrown. in cases when the current image does not
360    /// have the required tags on its own.
361    ///
362    /// ## Returns
363    /// The resolution as a lon-lat vector
364    pub fn resolution_ll(&self) -> VectorPoint<()> {
365        let mut resolution = self.resolution();
366        self.transformer.borrow_mut().forward_mut(&mut resolution);
367        resolution
368    }
369
370    /// Returns whether or not the pixels of the image depict an area (or point).
371    ///
372    /// ## Returns
373    /// Whether the pixels are a point
374    pub fn pixel_is_area(&self) -> bool {
375        self.image_directory
376            .geo_key_directory
377            .get_short(GeoKeyDirectoryKeys::GTRasterTypeGeoKey as u16)
378            .unwrap_or(0)
379            == 1
380    }
381
382    /// Returns the image bounding box as an array of 4 values: min-x, min-y,
383    /// max-x and max-y. When the image has no affine transformation, then an
384    /// exception is thrown.
385    ///
386    /// ## Parameters
387    /// - `transform`: apply affine transformation or proj4 transformation
388    ///
389    /// ## Returns
390    /// The bounding box
391    pub fn get_bbox(&mut self, transform: bool) -> BBox {
392        let height = self.height() as f64;
393        let width = self.width() as f64;
394        let model_transformation =
395            self.image_directory.variables.getf64s(FieldTagNames::ModelTransformation as u16);
396
397        if transform && let Some(model_transformation) = model_transformation {
398            let [a, b, _c, d, e, f, _g, h] = model_transformation[..8].try_into().unwrap();
399            let corners = [[0., 0.], [0., height], [width, 0.], [width, height]];
400            let projected = corners.map(|[_i, _j]| [d + a * _i + b * _j, h + e * _i + f * _j]);
401            let xs = projected.map(|pt| pt[0]);
402            let ys = projected.map(|pt| pt[1]);
403
404            BBox::new(
405                xs.iter().copied().fold(f64::INFINITY, fmin),
406                ys.iter().copied().fold(f64::INFINITY, fmin),
407                xs.iter().copied().fold(f64::NEG_INFINITY, fmax),
408                ys.iter().copied().fold(f64::NEG_INFINITY, fmax),
409            )
410        } else {
411            let VectorPoint { x: x1, y: y1, .. } = self.origin();
412            let VectorPoint { x: r1, y: r2, .. } = self.resolution();
413            let x2 = x1 + r1 * width;
414            let y2 = y1 + r2 * height;
415            let min_x = fmin(x1, x2);
416            let min_y = fmin(y1, y2);
417            let max_x = fmax(x1, x2);
418            let max_y = fmax(y1, y2);
419
420            if transform {
421                let mut min_vec: VectorPoint<()> = VectorPoint::new_xy(min_x, min_y, None);
422                self.transformer.borrow_mut().forward_mut(&mut min_vec);
423                let mut max_vec: VectorPoint<()> = VectorPoint::new_xy(max_x, max_y, None);
424                self.transformer.borrow_mut().forward_mut(&mut max_vec);
425                BBox::new(min_vec.x, min_vec.y, max_vec.x, max_vec.y)
426            } else {
427                BBox::new(min_x, min_y, max_x, max_y)
428            }
429        }
430    }
431
432    /// Returns the raster data of the image.
433    ///
434    /// ## Parameters
435    /// - `samples`: Samples to read from the image
436    ///
437    /// ## Returns
438    /// The raster data
439    pub fn raster_data(&mut self, samples: Option<Vec<u16>>) -> Raster {
440        let samples =
441            samples.unwrap_or_else(|| (0..self.samples_per_pixel()).map(|v| v as u16).collect());
442        // setup
443        let tile_width = self.tile_width();
444        let tile_height = self.tile_height();
445        let width = self.width();
446        let height = self.height();
447        let sample_per_pixel = self.samples_per_pixel();
448        let bits_per_sample = self
449            .image_directory
450            .variables
451            .get_u16s(FieldTagNames::BitsPerSample as u16)
452            .unwrap_or_default();
453        let mut bytes_per_pixel = self.bytes_per_pixel();
454        let SampleFormat { src_sample_offsets, sample_readers } =
455            self.get_sample_offsets_and_readers(&bits_per_sample, &samples);
456
457        let mut res: Vec<f64> = vec![0.0; width * height * sample_per_pixel];
458        let mut min = f64::INFINITY;
459        let mut max = f64::NEG_INFINITY;
460        let max_x_tile = width.div_ceil(tile_width);
461        let max_y_tile = height.div_ceil(tile_height);
462        for y_tile in 0..max_y_tile {
463            for x_tile in 0..max_x_tile {
464                let mut data: Option<Vec<u8>> = None;
465                if self.planar_configuration == 1 {
466                    data = Some(self.get_tile_or_strip(x_tile, y_tile, 0));
467                }
468                for sample_index in 0..samples.len() {
469                    let si = sample_index;
470                    let sample = samples[sample_index] as usize;
471                    if self.planar_configuration == 2 {
472                        bytes_per_pixel = bits_per_sample[sample].div_ceil(8) as usize;
473                        data = Some(self.get_tile_or_strip(x_tile, y_tile, sample));
474                    }
475                    let data = data.as_ref().expect("data failed to load");
476                    let block_height = self.block_height(y_tile);
477                    let first_line = y_tile * tile_height;
478                    let first_col = x_tile * tile_width;
479                    let last_line = first_line + block_height;
480                    let last_col = (x_tile + 1) * tile_width;
481                    let reader = sample_readers[si];
482
483                    let ymax = block_height
484                        .min(
485                            (block_height as isize - (last_line as isize - height as isize))
486                                as usize,
487                        )
488                        .min(height - first_line);
489                    let xmax = tile_width
490                        .min((tile_width as isize - (last_col as isize - width as isize)) as usize)
491                        .min(width - first_col);
492
493                    for y in 0..ymax {
494                        for x in 0..xmax {
495                            let pixel_offset = (y * tile_width + x) * bytes_per_pixel;
496                            let value = reader(
497                                data,
498                                pixel_offset + src_sample_offsets[si] as usize,
499                                self.little_endian,
500                            );
501                            if !value.is_nan() {
502                                min = min.min(value);
503                                max = max.max(value);
504                            }
505                            let window_coordinate = (y + first_line) * width * samples.len()
506                                + (x + first_col) * samples.len()
507                                + si;
508                            res[window_coordinate] = value;
509                        }
510                    }
511                }
512            }
513        }
514
515        Raster {
516            r#type: self.raster_array_type(),
517            data: res,
518            width,
519            height,
520            alpha: false,
521            min,
522            max,
523        }
524    }
525
526    /// Returns the RGBA raster data of the image.
527    ///
528    /// ## Returns
529    /// The RGBA raster data
530    pub fn get_rgba(&mut self) -> Raster {
531        let bits_per_sample =
532            self.image_directory.variables.get_u16s(FieldTagNames::BitsPerSample as u16);
533        let extra_samples = self
534            .image_directory
535            .variables
536            .get_u16s(FieldTagNames::ExtraSamples as u16)
537            .unwrap_or_else(|| vec![0])[0];
538        let pi: PhotometricInterpretations = self
539            .image_directory
540            .variables
541            .get_short(FieldTagNames::PhotometricInterpretation as u16)
542            .unwrap_or(0)
543            .into();
544        let samples = build_samples(pi, bits_per_sample, Some(extra_samples.into()));
545
546        let mut raster_data = self.raster_data(Some(samples));
547        let max = pow(2., self.get_bits_per_sample(None) as f64);
548        convert_color_space(
549            pi,
550            &mut raster_data,
551            max,
552            self.image_directory.variables.get_u16s(FieldTagNames::ColorMap as u16),
553        );
554        raster_data.alpha = extra_samples != 0;
555
556        raster_data
557    }
558
559    /// Build a vector feature from the image
560    ///
561    /// ## Returns
562    /// The vector feature with rgba values incoded into the points
563    pub fn get_multi_point_vector(&mut self) -> GeoTIFFVectorFeature {
564        let Raster { width, height, alpha, data, .. } = self.get_rgba();
565        let bbox = self.get_bbox(false);
566        let BBox { left: min_x, bottom: min_y, right: max_x, top: max_y } = bbox;
567        let mut coordinates: VectorMultiPoint<RGBA> = vec![];
568        let rgba_stride = if alpha { 4 } else { 3 };
569        let bound_x = if min_x == max_x { 1. } else { max_x - min_x };
570        let bound_y = if min_y == max_y { 1. } else { max_y - min_y };
571        let area_x_stride =
572            if self.pixel_is_area() { (0.5 / (width as f64)) * bound_x } else { 0. };
573        let area_y_stride =
574            if self.pixel_is_area() { (0.5 / (height as f64)) * bound_y } else { 0. };
575        let pixel_x_stride = if width == 1 { 1 } else { width - 1 };
576        let pixel_y_stride = if height == 1 { 1 } else { height - 1 };
577
578        for y in 0..height {
579            for x in 0..width {
580                // Adjust x_pos and y_pos relative to the bounding box
581                let x_pos =
582                    min_x + ((x as f64) / (pixel_x_stride as f64)) * bound_x + area_x_stride;
583                let y_pos =
584                    min_y + ((y as f64) / (pixel_y_stride as f64)) * bound_y + area_y_stride;
585                // Extract RGBA values
586                let r = data[y * width * rgba_stride + x * rgba_stride];
587                let g = data[y * width * rgba_stride + x * rgba_stride + 1];
588                let b = data[y * width * rgba_stride + x * rgba_stride + 2];
589                let a =
590                    if alpha { data[y * width * rgba_stride + x * rgba_stride + 3] } else { 255. };
591                // find the lon-lat coordinates of the point
592                let mut point: VectorPoint<RGBA> = VectorPoint::new_xy(
593                    x_pos,
594                    y_pos,
595                    Some(RGBA { r: r / 255., g: g / 255., b: b / 255., a: a / 255. }),
596                );
597                self.transformer.borrow_mut().forward_mut(&mut point);
598                // Add point to coordinates array
599                coordinates.push(point);
600            }
601        }
602
603        GeoTIFFVectorFeature::new_wm(
604            None,
605            Properties::default(),
606            VectorGeometry::new_multipoint(coordinates, Some(bbox.into())),
607            Some(GeoTIFFMetadata { width, height, alpha }),
608        )
609    }
610
611    /// Get the data for a tile or strip
612    ///
613    /// ## Parameters
614    /// - `x`: the tile or strip x coordinate
615    /// - `y`: the tile or strip y coordinate
616    /// - `sample`: the sample
617    ///
618    /// ## Returns
619    /// The data as a buffer
620    fn get_tile_or_strip(&mut self, x: usize, y: usize, sample: usize) -> Vec<u8> {
621        let num_tiles_per_row = self.width().div_ceil(self.tile_width());
622        let num_tiles_per_col = self.height().div_ceil(self.tile_height());
623        let index = if self.planar_configuration == 1 {
624            y * num_tiles_per_row + x
625        } else if self.planar_configuration == 2 {
626            sample * num_tiles_per_row * num_tiles_per_col + y * num_tiles_per_row + x
627        } else {
628            0
629        };
630
631        let tile_offsets = self
632            .image_directory
633            .variables
634            .get_u32s(FieldTagNames::TileOffsets as u16)
635            .unwrap_or_default();
636        let tile_byte_counts = self
637            .image_directory
638            .variables
639            .get_u32s(FieldTagNames::TileByteCounts as u16)
640            .unwrap_or_default();
641        let strip_offsets = self
642            .image_directory
643            .variables
644            .get_u32s(FieldTagNames::StripOffsets as u16)
645            .unwrap_or_default();
646        let strip_byte_counts = self
647            .image_directory
648            .variables
649            .get_u32s(FieldTagNames::StripByteCounts as u16)
650            .unwrap_or_default();
651
652        let offset = if self.is_tiled { tile_offsets.get(index) } else { strip_offsets.get(index) };
653        let byte_count =
654            if self.is_tiled { tile_byte_counts.get(index) } else { strip_byte_counts.get(index) };
655        if offset.is_none() || byte_count.is_none() {
656            panic!("Invalid offset or byte count");
657        }
658        let offset = *offset.unwrap() as u64;
659        let byte_count = *byte_count.unwrap() as u64;
660        let slice = self.reader.borrow_mut().slice(Some(offset), Some(offset + byte_count));
661        let mut data = if let Some(decode_fn) = &mut self.decode_fn {
662            decode_fn(
663                &slice,
664                self.image_directory
665                    .variables
666                    .get(FieldTagNames::JPEGTables as u16)
667                    .map(|v| v.0)
668                    .as_deref(),
669            )
670        } else {
671            slice
672        };
673        data = self.maybe_apply_predictor(data);
674        let sample_format = self.get_sample_format(None) as usize;
675        let bits_per_sample = self.get_bits_per_sample(None) as usize;
676
677        if needs_normalization(sample_format, bits_per_sample) {
678            normalize_array(
679                data,
680                sample_format,
681                self.planar_configuration as usize,
682                self.samples_per_pixel(),
683                bits_per_sample,
684                self.tile_width(),
685                self.block_height(y),
686            )
687        } else {
688            // convert data into f64
689            data
690        }
691    }
692
693    /// Apply the predictor if necessary
694    ///
695    /// ## Parameters
696    /// - `data`: the raw data
697    ///
698    /// ## Returns
699    /// The data with the predictor applied
700    fn maybe_apply_predictor(&mut self, data: Vec<u8>) -> Vec<u8> {
701        let predictor =
702            self.image_directory.variables.get_short(FieldTagNames::Predictor as u16).unwrap_or(1);
703        if predictor == 1 {
704            data
705        } else {
706            let tile_width = if self.is_tiled { self.tile_width() } else { self.width() };
707            let tile_height = if self.is_tiled {
708                self.tile_height()
709            } else {
710                let val = self
711                    .image_directory
712                    .variables
713                    .get_short(FieldTagNames::RowsPerStrip as u16)
714                    .unwrap_or_else(|| self.height() as i16);
715                val as usize
716            };
717            let bits_per_sample = self
718                .image_directory
719                .variables
720                .get_u16s(FieldTagNames::BitsPerSample as u16)
721                .unwrap_or_default();
722            apply_predictor(
723                data,
724                predictor,
725                tile_width,
726                tile_height,
727                bits_per_sample,
728                self.planar_configuration,
729            )
730        }
731    }
732
733    /// Get the sample format. If first time than build it
734    ///
735    /// ## Parameters
736    /// - `bits_per_sample`: the bits per sample
737    /// - `samples`: the samples
738    ///
739    /// ## Returns
740    /// The sample format
741    fn get_sample_offsets_and_readers(
742        &mut self,
743        bits_per_sample: &[u16],
744        samples: &[u16],
745    ) -> SampleFormat {
746        if self.src_sample_offsets.is_empty()
747            && self.sample_readers.is_empty()
748            && samples.len() == self.sample_readers.len()
749        {
750            return SampleFormat {
751                src_sample_offsets: self.src_sample_offsets.clone(),
752                sample_readers: self.sample_readers.clone(),
753            };
754        }
755        let mut src_sample_offsets: Vec<u32> = vec![];
756        let mut sample_readers = vec![];
757        for sample in samples.iter() {
758            if self.planar_configuration == 1 {
759                src_sample_offsets
760                    .push((sample_sum(bits_per_sample, 0, *sample as usize) / 8) as u32);
761            } else {
762                src_sample_offsets.push(0);
763            }
764            let format = self.get_sample_format(Some(*sample as usize));
765            sample_readers.push(get_reader_for_sample(*sample, format));
766        }
767        self.src_sample_offsets = src_sample_offsets.clone();
768        self.sample_readers = sample_readers.clone();
769
770        SampleFormat { src_sample_offsets, sample_readers }
771    }
772}
773
774//  /// Get a value in the image
775//  ///
776//  /// ## Parameters
777//  /// - `x`: the x coordinate
778//  /// - `y`: the y coordinate
779//  /// - `inv_y`: if true, the y coordinate is inverted
780//  /// - `samples`: Samples to read from the image
781//  ///
782//  /// ## Returns
783//  /// The sample
784//   pub fn get_value(
785//     x: number,
786//     y: number,
787//     inv_y = false,
788//     samples: number[] = [...Array(self.samples_per_pixel()).keys()],
789//   ): Promise<number[]> {
790//     // setup
791//     let { tile_width, tile_height, width } = this;
792//     let bits_per_sample = self.#imageDirectory.BitsPerSample ?? [];
793//     let bytesPerPixel = self.bytesPerPixel;
794//     let { src_sample_offsets, sample_readers } = self.get_sample_offsets_and_readers(
795//       bits_per_sample,
796//       samples,
797//     );
798//     let res: number[] = new Array(samples.length);
799//     // invert Y if needed
800//     if (inv_y) y = self.height - 1 - y;
801
802//     // search the right tile
803//     let x_tile = floor(x / tile_width);
804//     let y_tile = floor(y / tile_height);
805//     let data: ArrayBufferLike | undefined;
806//     if (self.planar_configuration === 1) {
807//       data = await self.get_tile_or_strip(x_tile, y_tile, 0);
808//     }
809//     for (let sample_index = 0; sample_index < samples.len(); ++sample_index) {
810//       let si = sample_index;
811//       let sample = samples[sample_index];
812//       if (self.planar_configuration === 2) {
813//         bytesPerPixel = Math.ceil(bits_per_sample[sample] / 8);
814//         data = await self.get_tile_or_strip(x_tile, y_tile, sample);
815//       }
816//       if (data === undefined) throw new Error('data failed to load');
817//       let data_view = new DataView(data);
818//       let first_line = y_tile * tile_height;
819//       let first_col = x_tile * tile_width;
820//       let reader = sample_readers[si];
821
822//       let pixel_offset = (y * tile_width + x) * bytesPerPixel;
823//       let value = reader.call(data_view, pixel_offset + src_sample_offsets[si], self.#little_endian);
824//       let window_coordinate =
825//         (y + first_line) * width * samples.len() + (x + first_col) * samples.len() + si;
826//       res[floor(window_coordinate) % samples.len()] = value;
827//     }
828
829//     return res;
830//   }
gistools/readers/geotiff/image.rs

gistools/readers/geotiff/
image.rs
