DParser Finished, should do more optimimzations

down.py

@@ -0,0 +1,197 @@
+def _downsample_2d(src, mask, use_mask, method, fill_value, mode_rank, out):
+    src_w = src.shape[-1]
+    src_h = src.shape[-2]
+    out_w = out.shape[-1]
+    out_h = out.shape[-2]
+
+    if src_w == out_w and src_h == out_h:
+        return src
+
+    if out_w > src_w or out_h > src_h:
+        raise ValueError("invalid target size")
+
+    scale_x = src_w / out_w
+    scale_y = src_h / out_h
+
+    if method == DS_FIRST or method == DS_LAST:
+        for out_y in range(out_h):
+            src_yf0 = scale_y * out_y
+            src_yf1 = src_yf0 + scale_y
+            src_y0 = int(src_yf0)
+            src_y1 = int(src_yf1)
+            if src_y1 == src_yf1 and src_y1 > src_y0:
+                src_y1 -= 1
+            for out_x in range(out_w):
+                src_xf0 = scale_x * out_x
+                src_xf1 = src_xf0 + scale_x
+                src_x0 = int(src_xf0)
+                src_x1 = int(src_xf1)
+                if src_x1 == src_xf1 and src_x1 > src_x0:
+                    src_x1 -= 1
+                done = False
+                value = fill_value
+                for src_y in range(src_y0, src_y1 + 1):
+                    for src_x in range(src_x0, src_x1 + 1):
+                        v = src[src_y, src_x]
+                        if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
+                            value = v
+                            if method == DS_FIRST:
+                                done = True
+                                break
+                    if done:
+                        break
+                out[out_y, out_x] = value
+
+    elif method == DS_MODE:
+        max_value_count = int(scale_x + 1) * int(scale_y + 1)
+        values = np.zeros((max_value_count,), dtype=src.dtype)
+        frequencies = np.zeros((max_value_count,), dtype=np.uint32)
+        for out_y in range(out_h):
+            src_yf0 = scale_y * out_y
+            src_yf1 = src_yf0 + scale_y
+            src_y0 = int(src_yf0)
+            src_y1 = int(src_yf1)
+            wy0 = 1.0 - (src_yf0 - src_y0)
+            wy1 = src_yf1 - src_y1
+            if wy1 < _EPS:
+                wy1 = 1.0
+                if src_y1 > src_y0:
+                    src_y1 -= 1
+            for out_x in range(out_w):
+                src_xf0 = scale_x * out_x
+                src_xf1 = src_xf0 + scale_x
+                src_x0 = int(src_xf0)
+                src_x1 = int(src_xf1)
+                wx0 = 1.0 - (src_xf0 - src_x0)
+                wx1 = src_xf1 - src_x1
+                if wx1 < _EPS:
+                    wx1 = 1.0
+                    if src_x1 > src_x0:
+                        src_x1 -= 1
+                value_count = 0
+                for src_y in range(src_y0, src_y1 + 1):
+                    wy = wy0 if (src_y == src_y0) else wy1 if (src_y == src_y1) else 1.0
+                    for src_x in range(src_x0, src_x1 + 1):
+                        wx = wx0 if (src_x == src_x0) else wx1 if (src_x == src_x1) else 1.0
+                        v = src[src_y, src_x]
+                        if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
+                            w = wx * wy
+                            found = False
+                            for i in range(value_count):
+                                if v == values[i]:
+                                    frequencies[i] += w
+                                    found = True
+                                    break
+                            if not found:
+                                values[value_count] = v
+                                frequencies[value_count] = w
+                                value_count += 1
+                w_max = -1.
+                value = fill_value
+                if mode_rank == 1:
+                    for i in range(value_count):
+                        w = frequencies[i]
+                        if w > w_max:
+                            w_max = w
+                            value = values[i]
+                elif mode_rank <= max_value_count:
+                    max_frequencies = np.full(mode_rank, -1.0, dtype=np.float64)
+                    indices = np.zeros(mode_rank, dtype=np.int64)
+                    for i in range(value_count):
+                        w = frequencies[i]
+                        for j in range(mode_rank):
+                            if w > max_frequencies[j]:
+                                max_frequencies[j] = w
+                                indices[j] = i
+                                break
+                    value = values[indices[mode_rank - 1]]
+
+                out[out_y, out_x] = value
+
+    elif method == DS_MEAN:
+        for out_y in range(out_h):
+            src_yf0 = scale_y * out_y
+            src_yf1 = src_yf0 + scale_y
+            src_y0 = int(src_yf0)
+            src_y1 = int(src_yf1)
+            wy0 = 1.0 - (src_yf0 - src_y0)
+            wy1 = src_yf1 - src_y1
+            if wy1 < _EPS:
+                wy1 = 1.0
+                if src_y1 > src_y0:
+                    src_y1 -= 1
+            for out_x in range(out_w):
+                src_xf0 = scale_x * out_x
+                src_xf1 = src_xf0 + scale_x
+                src_x0 = int(src_xf0)
+                src_x1 = int(src_xf1)
+                wx0 = 1.0 - (src_xf0 - src_x0)
+                wx1 = src_xf1 - src_x1
+                if wx1 < _EPS:
+                    wx1 = 1.0
+                    if src_x1 > src_x0:
+                        src_x1 -= 1
+                v_sum = 0.0
+                w_sum = 0.0
+                for src_y in range(src_y0, src_y1 + 1):
+                    wy = wy0 if (src_y == src_y0) else wy1 if (src_y == src_y1) else 1.0
+                    for src_x in range(src_x0, src_x1 + 1):
+                        wx = wx0 if (src_x == src_x0) else wx1 if (src_x == src_x1) else 1.0
+                        v = src[src_y, src_x]
+                        if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
+                            w = wx * wy
+                            v_sum += w * v
+                            w_sum += w
+                if w_sum < _EPS:
+                    out[out_y, out_x] = fill_value
+                else:
+                    out[out_y, out_x] = v_sum / w_sum
+
+    elif method == DS_VAR or method == DS_STD:
+        for out_y in range(out_h):
+            src_yf0 = scale_y * out_y
+            src_yf1 = src_yf0 + scale_y
+            src_y0 = int(src_yf0)
+            src_y1 = int(src_yf1)
+            wy0 = 1.0 - (src_yf0 - src_y0)
+            wy1 = src_yf1 - src_y1
+            if wy1 < _EPS:
+                wy1 = 1.0
+                if src_y1 > src_y0:
+                    src_y1 -= 1
+            for out_x in range(out_w):
+                src_xf0 = scale_x * out_x
+                src_xf1 = src_xf0 + scale_x
+                src_x0 = int(src_xf0)
+                src_x1 = int(src_xf1)
+                wx0 = 1.0 - (src_xf0 - src_x0)
+                wx1 = src_xf1 - src_x1
+                if wx1 < _EPS:
+                    wx1 = 1.0
+                    if src_x1 > src_x0:
+                        src_x1 -= 1
+                v_sum = 0.0
+                w_sum = 0.0
+                wv_sum = 0.0
+                wvv_sum = 0.0
+                for src_y in range(src_y0, src_y1 + 1):
+                    wy = wy0 if (src_y == src_y0) else wy1 if (src_y == src_y1) else 1.0
+                    for src_x in range(src_x0, src_x1 + 1):
+                        wx = wx0 if (src_x == src_x0) else wx1 if (src_x == src_x1) else 1.0
+                        v = src[src_y, src_x]
+                        if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
+                            w = wx * wy
+                            v_sum += v
+                            w_sum += w
+                            wv_sum += w * v
+                            wvv_sum += w * v * v
+                if w_sum < _EPS:
+                    out[out_y, out_x] = fill_value
+                else:
+                    out[out_y, out_x] = (wvv_sum * w_sum - wv_sum * wv_sum) / w_sum / w_sum
+        if method == DS_STD:
+            out = np.sqrt(out).astype(out.dtype)
+    else:
+        raise ValueError('invalid upsampling method')
+
+    return out

src/data/mapper.rs

@@ -13,6 +13,7 @@ pub struct Mapper {
     range: (Range<f64>, Range<f64>),
     bounds: (f64, f64, f64, f64),
 }
+unsafe impl Sync for Mapper {}
 
 impl From<Proj> for Mapper {
     fn from(proj: Proj) -> Self {
@@ -77,15 +78,15 @@ impl Mapper {
     pub fn ring_map(&self, ring: &LineString) -> Result<LineString, ProjError> {
         let mut result = Vec::new();
         for l in ring.lines() {
-            let start_projected = self.map((l.start.x, l.start.y))?;
+            let start_projected: (f64, f64) = l.start.into();
-            let end_projected = self.map((l.end.x, l.end.y))?;
+            let end_projected: (f64, f64) = l.end.into();
-
             let cartesian_start = self.cartesian(l.start.x, l.start.y);
             let cartesian_end = self.cartesian(l.end.x, l.end.y);
 
             let delta2 = 0.5;
             let depth = 16;
             let mut res: Vec<GCoord> = Vec::new();
+            res.push(l.start);
             self.resample_line_to(
                 start_projected,
                 end_projected,
@@ -98,6 +99,7 @@ impl Mapper {
                 depth,
                 &mut res,
             )?;
+            res.push(l.end);
             result.extend(res);
         }
 

src/data/mod.rs

@@ -1,5 +1,5 @@
 use ndarray::{
-    s, Array, Array1, Array2, Array3, ArrayBase, Ix1, Ix2, OwnedRepr, RawDataClone, ViewRepr,
+    s, Array, Array1, Array2, Array3, ArrayBase, Axis, Ix1, Ix2, OwnedRepr, RawDataClone, ViewRepr,
 };
 use npyz::{npz::NpzArchive, Deserialize};
 use num_traits::{AsPrimitive, FromPrimitive, Num};
@@ -21,11 +21,10 @@ pub enum DataError {
     },
 }
 
-#[derive(PartialEq)]
+pub enum DownSampleMeth {
-pub enum Axis {
+    STD,
-    Asix0,
+    MEAN,
-    Axis1,
+    VAR,
-    Axis2,
 }
 
 #[derive(Clone, Copy)]
@@ -48,6 +47,16 @@ where
     pub coord_type: CoordType,
 }
 
+pub type Radar2d<T> = RadarData2d<T, OwnedRepr<T>>;
+
+impl<T: Num + Clone + PartialEq + PartialOrd> Radar2d<T> {
+    pub fn load(path: impl AsRef<Path>, meth: impl DataLoader<Self>) -> Result<Self, DataError> {
+        Ok(meth.load(path)?)
+    }
+}
+
+pub type Radar2dRef<'a, T> = RadarData2d<T, ViewRepr<&'a T>>;
+
 pub struct RadarData3d<T, X = f64, Y = f64, Z = f64>
 where
     T: Num,
@@ -68,58 +77,134 @@ where
     T: Num + AsPrimitive<f64> + FromPrimitive + Clone + Debug + PartialOrd + PartialEq,
     Raw: ndarray::Data<Elem = T> + Clone + ndarray::RawDataClone,
 {
-    fn value_to_owned(self) -> RadarData2d<T, OwnedRepr<T>> {
+    pub fn downsample(&self, output_shape: (usize, usize), meth: DownSampleMeth) -> Radar2d<T> {
-        RadarData2d {
+        let shape = self.data.shape();
-            dim1: self.dim1,
+        assert!(shape[0] > output_shape.0 && shape[1] > output_shape.1);
-            dim2: self.dim2,
+
+        if shape[0] == output_shape.0 && shape[1] == output_shape.1 {
+            return Radar2d {
+                dim1: self.dim1.to_owned(),
+                dim2: self.dim2.to_owned(),
                 data: self.data.to_owned(),
-            coord_type: self.coord_type,
+                coord_type: self.coord_type.clone(),
-        }
+            };
         }
 
-    fn resample(
+        let scale_x = shape[1] as f64 / output_shape.1 as f64;
+        let scale_y = shape[0] as f64 / output_shape.0 as f64;
+
+        let mut output: Array2<T> = Array2::zeros(output_shape);
+
+        let mut dim1 = Array1::zeros(output_shape.1);
+        let mut dim2 = Array1::zeros(output_shape.0);
+
+        output.iter_mut().enumerate().for_each(|(s, vv)| {
+            let ri = s / output_shape.1;
+            let ci = s % output_shape.1;
+            let src_yf0 = scale_y * ci as f64;
+            let src_yf1 = src_yf0 + scale_y;
+            let src_y0 = src_yf0.floor() as usize;
+            let src_y1 = src_yf1.floor() as usize;
+            let src_xf0 = scale_x * ri as f64;
+            let src_xf1 = src_xf0 + scale_x;
+            let src_x0 = src_xf0.floor() as usize;
+            let src_x1 = src_xf1.floor() as usize;
+
+            match meth {
+                DownSampleMeth::MEAN => {}
+                DownSampleMeth::VAR | DownSampleMeth::STD => self.var(
+                    src_yf0, src_yf1, src_xf0, src_xf1, src_y0, src_y1, src_x0, src_x1, vv,
+                ),
+            }
+
+            dim1[ci] = self.dim1[src_x0];
+            dim2[ri] = self.dim2[src_y0];
+        });
+
+        if let DownSampleMeth::STD = meth {
+            output = output.mapv(|x| T::from_f64(f64::sqrt(x.as_())).unwrap());
+        }
+
+        return Radar2d {
+            dim1,
+            dim2,
+            data: output,
+            coord_type: self.coord_type.clone(),
+        };
+    }
+
+    fn var(
         &self,
-        width_rate: f64,
+        src_yf0: f64,
-        height_rate: f64,
+        src_yf1: f64,
-        filter_len: f64,
+        src_xf0: f64,
-    ) -> Result<RadarData2d<T, OwnedRepr<T>>, DataError> {
+        src_xf1: f64,
-        let width_rate = width_rate.min(1.0);
+        src_y0: usize,
-        let height_rate = height_rate.min(1.0);
+        mut src_y1: usize,
-        match self.coord_type {
+        src_x0: usize,
-            CoordType::Polar => Err(DataError::FormatError),
+        mut src_x1: usize,
-            CoordType::LatLon => {
+        row: &mut T,
-                let width_filtered: ArrayBase<ndarray::OwnedRepr<T>, Ix2> =
+    ) {
-                    Self::_resample(&self.data, width_rate, filter_len);
+        let wy0 = 1f64 - (src_yf0 - src_y0 as f64);
+        let mut wy1 = src_yf1 - src_y1 as f64;
+        let wx0 = 1f64 - (src_xf0 - src_x0 as f64);
+        let mut wx1 = src_xf1 - src_x1 as f64;
 
-                let result: ArrayBase<OwnedRepr<T>, Ix2> =
+        if wy1 < f64::EPSILON {
-                    Self::_resample(&width_filtered.t(), height_rate, filter_len)
+            wy1 = 1f64;
-                        .t()
+            if src_y1 > src_y0 {
-                        .to_owned();
+                src_y1 -= 1;
-
-                let new_dim1: ArrayBase<OwnedRepr<f64>, Ix1> = Self::_resample(
-                    &Array2::from_shape_vec((1, self.dim1.len()), self.dim1.to_vec()).unwrap(),
-                    width_rate,
-                    filter_len,
-                )
-                .slice(s![0, ..])
-                .to_owned();
-
-                let new_dim2: ArrayBase<OwnedRepr<f64>, Ix1> = Self::_resample(
-                    &Array2::from_shape_vec((1, self.dim2.len()), self.dim2.to_vec()).unwrap(),
-                    height_rate,
-                    filter_len,
-                )
-                .slice(s![0, ..])
-                .to_owned();
-
-                Ok(RadarData2d {
-                    dim1: new_dim1,
-                    dim2: new_dim2,
-                    data: result,
-                    coord_type: self.coord_type.to_owned(),
-                })
             }
         }
+
+        if wx1 < f64::EPSILON {
+            wx1 = 1f64;
+            if src_x1 > src_x0 {
+                src_x1 -= 1;
+            }
+        }
+
+        let mut v_sum: f64 = 0.0;
+        let mut w_sum = 0f64;
+        let mut wv_sum = 0f64;
+        let mut wvv_sum = 0f64;
+
+        for src_y in src_y0..=src_y1 {
+            let wy = if src_y == src_y0 {
+                wy0
+            } else {
+                if src_y == src_y1 {
+                    wy1
+                } else {
+                    1.0
+                }
+            };
+            for src_x in src_x0..=src_x1 {
+                let wx = if src_x == src_x0 {
+                    wx0
+                } else {
+                    if src_x == src_x1 {
+                        wx1
+                    } else {
+                        1.0
+                    }
+                };
+
+                let v = self.data[[src_y, src_x]].as_();
+                let w = wx * wy;
+
+                v_sum += v;
+                w_sum += w;
+                wv_sum += w * v;
+                wvv_sum += w * v * v;
+            }
+
+            if w_sum < f64::EPSILON {
+                return;
+            }
+
+            *row = T::from_f64((wvv_sum * w_sum - wv_sum * wv_sum) / w_sum / w_sum).unwrap();
+        }
     }
 
     fn split(&self) -> [RadarData2d<T, ndarray::ViewRepr<&T>>; 4] {
@@ -158,121 +243,6 @@ where
             },
         ];
     }
-
-    fn _resample<'a, V, R: ndarray::Data<Elem = V>>(
-        data: &'a ArrayBase<R, Ix2>,
-        rate: f64,
-        filter_len: f64,
-    ) -> Array2<V>
-    where
-        V: Num + Clone + AsPrimitive<f64> + FromPrimitive,
-    {
-        let ori_width = data.ncols();
-        let ori_height = data.nrows();
-
-        let new_width = (ori_width as f64 * rate).ceil() as usize;
-        let mut result: Array2<V> = Array2::zeros((ori_height, new_width));
-        (0..ori_height).into_iter().for_each(|height| {
-            for width in 0..new_width {
-                let center_x = (width as f64 + 0.5) / new_width as f64 * ori_width as f64;
-                let filter_start = center_x - filter_len / 2.0;
-                let start_idx = (filter_start - 0.5).ceil() as usize;
-
-                let mut value_sum = 0.0;
-                let mut filter_sum = 0.0;
-
-                for i in 0..filter_len as usize {
-                    let input_x = start_idx + i;
-                    let weight = windowed_sinc(
-                        (input_x as f64 + 0.5 - center_x) * rate,
-                        (input_x as f64 + 0.5 - filter_start) / filter_len,
-                    );
-                    value_sum += weight * data[[height, input_x.clamp(0, ori_width - 1)]].as_();
-                    filter_sum += weight;
-                }
-
-                result[[height, width]] = V::from_f64(value_sum / filter_sum).unwrap();
-            }
-        });
-        result
-    }
-}
-
-#[inline]
-fn windowed_sinc(x: f64, y: f64) -> f64 {
-    let x = x * PI;
-    let sinc = if x != 0.0 { f64::sin(x) / x } else { 1.0 };
-    let window = if 0f64 <= y && y <= 1.0 {
-        1.0 - (y - 0.5).abs() * 2.0
-    } else {
-        0f64
-    };
-    sinc * window
-}
-
-pub struct LevelData<T: Num + Clone + PartialEq + PartialOrd>(
-    pub Quadtree<i64, RadarData2d<T, OwnedRepr<T>>>,
-);
-
-impl<T> LevelData<T>
-where
-    T: Num + Clone + AsPrimitive<f64> + FromPrimitive + Debug + PartialEq + PartialOrd,
-{
-    fn value(
-        level_data: Vec<RadarData2d<T, OwnedRepr<T>>>,
-        level_num: usize,
-    ) -> Vec<RadarData2d<T, OwnedRepr<T>>> {
-        if level_num == 0 {
-            return level_data;
-        }
-
-        let mut result: Vec<RadarData2d<T, OwnedRepr<T>>> =
-            Vec::with_capacity(level_data.len() * 4);
-
-        let results = level_data
-            .iter()
-            .flat_map(|v| v.split().into_iter().map(|x| x.value_to_owned()));
-
-        result.extend(results);
-
-        return Self::value(result, level_num - 1);
-    }
-
-    fn new(data: &RadarData2d<T, OwnedRepr<T>>, level: usize, rate: f64) -> Self {
-        // let rate = 1.0 / level as f64;
-        let resampled = data.resample(rate, rate, 2.0).unwrap();
-        let blocks = Self::value(vec![resampled], level);
-        let mut tree: Quadtree<i64, RadarData2d<T, OwnedRepr<T>>> =
-            quadtree_rs::Quadtree::new(level);
-
-        blocks.into_iter().for_each(|block| {
-            tree.insert(
-                AreaBuilder::default()
-                    .anchor(quadtree_rs::point::Point {
-                        x: *block.dim1.first().unwrap() as i64,
-                        y: *block.dim2.first().unwrap() as i64,
-                    })
-                    .dimensions((block.dim1.len() as i64, block.dim2.len() as i64))
-                    .build()
-                    .unwrap(),
-                block,
-            );
-        });
-
-        Self(tree)
-    }
-}
-
-pub fn levels<T>(data: Radar2d<T>, levels: usize) -> Vec<LevelData<T>>
-where
-    T: Num + Clone + AsPrimitive<f64> + FromPrimitive + Debug,
-    T: PartialEq + PartialOrd,
-{
-    let numerator = 1.0 / levels as f64;
-    (0..levels)
-        .into_iter()
-        .map(|level| LevelData::new(&data, level + 1, 1.0 - level as f64 * numerator))
-        .collect()
 }
 
 impl<T, Raw> MultiDimensionData for RadarData2d<T, Raw>
@@ -340,12 +310,162 @@ where
     }
 }
 
-pub type Radar2d<T> = RadarData2d<T, OwnedRepr<T>>;
+// fn resample(
+//     &self,
+//     width_rate: f64,
+//     height_rate: f64,
+//     filter_len: f64,
+// ) -> Result<RadarData2d<T, OwnedRepr<T>>, DataError> {
+//     let width_rate = width_rate.min(1.0);
+//     let height_rate = height_rate.min(1.0);
+//     match self.coord_type {
+//         CoordType::Polar => Err(DataError::FormatError),
+//         CoordType::LatLon => {
+//             let width_filtered: ArrayBase<ndarray::OwnedRepr<T>, Ix2> =
+//                 Self::_resample(&self.data, width_rate, filter_len);
 
-impl<T: Num + Clone + PartialEq + PartialOrd> Radar2d<T> {
+//             let result: ArrayBase<OwnedRepr<T>, Ix2> =
-    pub fn load(path: impl AsRef<Path>, meth: impl DataLoader<Self>) -> Result<Self, DataError> {
+//                 Self::_resample(&width_filtered.t(), height_rate, filter_len)
-        Ok(meth.load(path)?)
+//                     .t()
-    }
+//                     .to_owned();
+
+//             let new_dim1: ArrayBase<OwnedRepr<f64>, Ix1> = Self::_resample(
+//                 &Array2::from_shape_vec((1, self.dim1.len()), self.dim1.to_vec()).unwrap(),
+//                 width_rate,
+//                 filter_len,
+//             )
+//             .slice(s![0, ..])
+//             .to_owned();
+
+//             let new_dim2: ArrayBase<OwnedRepr<f64>, Ix1> = Self::_resample(
+//                 &Array2::from_shape_vec((1, self.dim2.len()), self.dim2.to_vec()).unwrap(),
+//                 height_rate,
+//                 filter_len,
+//             )
+//             .slice(s![0, ..])
+//             .to_owned();
+
+//             Ok(RadarData2d {
+//                 dim1: new_dim1,
+//                 dim2: new_dim2,
+//                 data: result,
+//                 coord_type: self.coord_type.to_owned(),
+//             })
+//         }
+//     }
+// }
+
+// fn _resample<'a, V, R: ndarray::Data<Elem = V>>(
+//         data: &'a ArrayBase<R, Ix2>,
+//         rate: f64,
+//         filter_len: f64,
+//     ) -> Array2<V>
+//     where
+//         V: Num + Clone + AsPrimitive<f64> + FromPrimitive,
+//     {
+//         let ori_width = data.ncols();
+//         let ori_height = data.nrows();
+
+//         let new_width = (ori_width as f64 * rate).ceil() as usize;
+//         let mut result: Array2<V> = Array2::zeros((ori_height, new_width));
+//         (0..ori_height).into_iter().for_each(|height| {
+//             for width in 0..new_width {
+//                 let center_x = (width as f64 + 0.5) / new_width as f64 * ori_width as f64;
+//                 let filter_start = center_x - filter_len / 2.0;
+//                 let start_idx = (filter_start - 0.5).ceil() as usize;
+
+//                 let mut value_sum = 0.0;
+//                 let mut filter_sum = 0.0;
+
+//                 for i in 0..filter_len as usize {
+//                     let input_x = start_idx + i;
+//                     let weight = windowed_sinc(
+//                         (input_x as f64 + 0.5 - center_x) * rate,
+//                         (input_x as f64 + 0.5 - filter_start) / filter_len,
+//                     );
+//                     value_sum += weight * data[[height, input_x.clamp(0, ori_width - 1)]].as_();
+//                     filter_sum += weight;
+//                 }
+
+//                 result[[height, width]] = V::from_f64(value_sum / filter_sum).unwrap();
+//             }
+//         });
+//         result
+//     }
+
+// pub struct LevelData<T: Num + Clone + PartialEq + PartialOrd>(
+//     pub Quadtree<i64, RadarData2d<T, OwnedRepr<T>>>,
+// );
+
+// impl<T> LevelData<T>
+// where
+//     T: Num + Clone + AsPrimitive<f64> + FromPrimitive + Debug + PartialEq + PartialOrd,
+// {
+//     fn value(
+//         level_data: Vec<RadarData2d<T, OwnedRepr<T>>>,
+//         level_num: usize,
+//     ) -> Vec<RadarData2d<T, OwnedRepr<T>>> {
+//         if level_num == 0 {
+//             return level_data;
+//         }
+
+//         let mut result: Vec<RadarData2d<T, OwnedRepr<T>>> =
+//             Vec::with_capacity(level_data.len() * 4);
+
+//         let results = level_data
+//             .iter()
+//             .flat_map(|v| v.split().into_iter().map(|x| x.value_to_owned()));
+
+//         result.extend(results);
+
+//         return Self::value(result, level_num - 1);
+//     }
+
+//     fn new(data: &RadarData2d<T, OwnedRepr<T>>, level: usize, rate: f64) -> Self {
+//         // let rate = 1.0 / level as f64;
+//         let resampled = data.resample(rate, rate, 2.0).unwrap();
+//         let blocks = Self::value(vec![resampled], level);
+//         let mut tree: Quadtree<i64, RadarData2d<T, OwnedRepr<T>>> =
+//             quadtree_rs::Quadtree::new(level);
+
+//         blocks.into_iter().for_each(|block| {
+//             tree.insert(
+//                 AreaBuilder::default()
+//                     .anchor(quadtree_rs::point::Point {
+//                         x: *block.dim1.first().unwrap() as i64,
+//                         y: *block.dim2.first().unwrap() as i64,
+//                     })
+//                     .dimensions((block.dim1.len() as i64, block.dim2.len() as i64))
+//                     .build()
+//                     .unwrap(),
+//                 block,
+//             );
+//         });
+
+//         Self(tree)
+//     }
+// }
+
+// pub fn levels<T>(data: Radar2d<T>, levels: usize) -> Vec<LevelData<T>>
+// where
+//     T: Num + Clone + AsPrimitive<f64> + FromPrimitive + Debug,
+//     T: PartialEq + PartialOrd,
+// {
+//     let numerator = 1.0 / levels as f64;
+//     (0..levels)
+//         .into_iter()
+//         .map(|level| LevelData::new(&data, level + 1, 1.0 - level as f64 * numerator))
+//         .collect()
+// }
+
+#[inline]
+fn windowed_sinc(x: f64, y: f64) -> f64 {
+    let x = x * PI;
+    let sinc = if x != 0.0 { f64::sin(x) / x } else { 1.0 };
+    let window = if 0f64 <= y && y <= 1.0 {
+        1.0 - (y - 0.5).abs() * 2.0
+    } else {
+        0f64
+    };
+    sinc * window
 }
-
-pub type Radar2dRef<T> = RadarData2d<T, ViewRepr<T>>;

src/main.rs

@@ -1,12 +1,7 @@
 use data::mapper::Mapper;
 use data::{Npz, Radar2d};
-use glib::{timeout_add, timeout_add_local};
+use gtk::prelude::*;
-use glib_macros::clone;
+use gtk::{gio, glib, Application, ApplicationWindow};
-use gtk::{
-    gio, glib, style_context_add_provider_for_display, Application, ApplicationWindow, CssProvider,
-    FontButton, GestureDrag,
-};
-use gtk::{prelude::*, DrawingArea};
 use std::ptr;
 mod data;
 mod monitor;
@@ -15,7 +10,6 @@ mod render;
 mod tree;
 mod window;
 use monitor::Monitor;
-use ndarray::parallel::prelude::{IntoParallelIterator, ParallelIterator};
 use painter::wgs84::{LatLonCoord, Mercator, ProjectionS, Range};
 use render::{BackgroundConfig, BackgroundWidget, ForegroundConfig, ForegroundWidget, Render};
 
@@ -63,17 +57,19 @@ fn build_ui(app: &Application) {
     let foreground_widget = ForegroundWidget::new(foreground_config);
     let render = Render::new(background_widget, foreground_widget);
 
-    let path = "/Users/ruomu/test2.npz";
+    let path = "/home/tsuki/projects/radar-g/test2.npz";
 
     let data = Radar2d::<i8>::load(path, Npz).unwrap();
 
     let projection = Mercator::new();
-    let mapper: Mapper = projection.into();
+    let mut mapper: Mapper = projection.into();
+    mapper.set_lat_range(29.960..30.764);
+    mapper.set_lon_range(120.038..120.965);
 
     render.set_mapper(mapper);
     let monitor = Monitor::new(render);
 
-    monitor.load_data_2d(data);
+    monitor.load_data_2d(data).unwrap();
 
     window.set_child(Some(&monitor));
     window.set_default_width(1000);

src/monitor/mod.rs

@@ -2,13 +2,14 @@ mod imp;
 use crate::data::RadarData2d;
 use crate::render::Render;
 use glib::clone;
-use glib::{clone::Downgrade, subclass::prelude::*};
+use glib::subclass::prelude::*;
 use gtk::glib;
 use gtk::traits::WidgetExt;
-use gtk::{EventController, EventControllerScrollFlags, Inhibit};
+use gtk::{EventControllerScrollFlags, Inhibit};
-use num_traits::Num;
+use num_traits::{AsPrimitive, FromPrimitive, Num};
 use proj::{Proj, ProjError};
 use std::borrow::Borrow;
+use std::fmt::Debug;
 
 glib::wrapper! {
     pub struct Monitor(ObjectSubclass<imp::Monitor>)
@@ -52,7 +53,14 @@ impl Monitor {
 
     pub fn load_data_2d<T, Raw>(&self, data: RadarData2d<T, Raw>) -> Result<(), ProjError>
     where
-        T: Num + Clone + PartialEq + PartialOrd,
+        T: Num
+            + Clone
+            + PartialEq
+            + PartialOrd
+            + AsPrimitive<i8>
+            + AsPrimitive<f64>
+            + Debug
+            + FromPrimitive,
         Raw: ndarray::Data<Elem = T> + Clone + ndarray::RawDataClone,
     {
         let renderer = self.imp().renderer.borrow();

src/render/foreground/imp.rs

@@ -29,6 +29,7 @@ pub struct ForegroundWidget {
     pub(super) config: RefCell<ForegroundConfig>,
     pub(super) dim1: RefCell<Option<Array2<f64>>>,
     pub(super) dim2: RefCell<Option<Array2<f64>>>,
+    pub(super) data: RefCell<Option<Array2<i8>>>,
 }
 
 #[glib::object_subclass]

src/render/foreground/mod.rs

@@ -1,9 +1,15 @@
 mod imp;
+use crate::tree::get;
 use crate::{data::mapper::Mapper, render::WindowCoord};
 use femtovg::{renderer::OpenGl, Canvas, Path};
+use femtovg::{Color, Paint};
+use geo_types::LineString;
 use glib::subclass::types::ObjectSubclassIsExt;
 use gtk::{ffi::gtk_widget_get_width, glib, graphene::Rect, prelude::SnapshotExtManual};
-use ndarray::Array2;
+use ndarray::parallel;
+use ndarray::Slice;
+use ndarray::{s, Array2, Axis, Zip};
+use std::cell::Ref;
 use std::ops::Range;
 
 pub use self::imp::ForegroundConfig;
@@ -26,14 +32,101 @@ impl ForegroundWidget {
         this
     }
 
-    pub fn draw(&self, canvas: &mut Canvas<OpenGl>, scale: f32, dpi: i32, mapper: &Mapper) {
+    pub fn draw(&self, canvas: &mut Canvas<OpenGl>, scale: f32, dpi: i32, mapper: Ref<'_, Mapper>) {
-        let canvas_widht = canvas.width();
+        let canvas_width = canvas.width();
         let canvas_height = canvas.height();
         let config = self.imp().config.borrow();
+
+        let dim1 = self.imp().dim1.borrow();
+        let dim2 = self.imp().dim2.borrow();
+        let data = self.imp().data.borrow();
+
+        let c = dim1.as_ref().unwrap().slice(s![..;3,..;3]);
+        let d = dim2.as_ref().unwrap().slice(s![..;3,..;3]);
+        let e = data.as_ref().unwrap().slice(s![..;3,..;3]);
+
+        let levels: Vec<i8> = vec![0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65];
+
+        let colors = vec![
+            Color::rgb(0, 172, 164),
+            Color::rgb(192, 192, 254),
+            Color::rgb(122, 114, 238),
+            Color::rgb(30, 38, 208),
+            Color::rgb(166, 252, 168),
+            Color::rgb(0, 234, 0),
+            Color::rgb(16, 146, 26),
+            Color::rgb(252, 244, 100),
+            Color::rgb(200, 200, 2),
+            Color::rgb(140, 140, 0),
+            Color::rgb(254, 172, 172),
+            Color::rgb(254, 100, 92),
+            Color::rgb(238, 2, 48),
+            Color::rgb(212, 142, 254),
+            Color::rgb(170, 36, 250),
+        ];
+
+        Zip::from(c).and(d).and(e).for_each(|lon, lat, d| {
+            if *d < -5 || *d > 70 {
+                return;
+            }
+
+            let left_bottom = mapper.map((*lon, *lat)).unwrap();
+            let right_top = mapper.map((lon + 0.003, lat + 0.003)).unwrap();
+
+            // let (lon1, lat1) = (
+            //     left_bottom.0 * canvas_width as f64,
+            //     left_bottom.1 * canvas_height as f64,
+            // );
+
+            let color = get(&levels, &colors, *d);
+            //
+            // let (lon2, lat2) = (
+            //     right_top.0 * canvas_width as f64,
+            //     right_top.1 * canvas_height as f64,
+            // );
+
+            let line_string: LineString = vec![
+                (left_bottom.0, left_bottom.1),
+                (right_top.0, left_bottom.1),
+                (right_top.0, right_top.1),
+                (left_bottom.0, right_top.1),
+                (left_bottom.0, left_bottom.1),
+            ]
+            .into();
+
+            // let line_string: LineString = vec![
+            //     (lon1, lat1),
+            //     (lon2, lat1),
+            //     (lon2, lat2),
+            //     (lon1, lat2),
+            //     (lon1, lat1),
+            // ]
+            // .into();
+
+            let res = mapper.ring_map(&line_string).unwrap();
+
+            let mut path = Path::new();
+
+            path.move_to(
+                left_bottom.0 as f32 * canvas_width,
+                left_bottom.1 as f32 * canvas_height,
+            );
+
+            for p in line_string.points() {
+                // path.move_to(p.x() as f32 * canvas_width, p.y() as f32 * canvas_height);
+                path.line_to(p.x() as f32 * canvas_width, p.y() as f32 * canvas_height);
+            }
+            canvas.fill_path(&path, &Paint::color(color));
+        });
+        canvas.flush();
     }
 
     pub(super) fn set_dims(&mut self, dims: (Array2<f64>, Array2<f64>)) {
         self.imp().dim1.replace(Some(dims.0));
         self.imp().dim2.replace(Some(dims.1));
     }
+
+    pub(super) fn set_data(&mut self, data: Array2<i8>) {
+        self.imp().data.replace(Some(data));
+    }
 }

src/render/imp.rs

@@ -51,10 +51,10 @@ impl ObjectImpl for Render {
         self.parent_constructed();
         let area = self.obj();
         area.set_has_stencil_buffer(true);
-        area.add_tick_callback(|area, _| {
+        // area.add_tick_callback(|area, _| {
-            area.queue_render();
+        //     area.queue_render();
-            glib::Continue(true)
+        //     glib::Continue(true)
-        });
+        // });
     }
 }
 
@@ -96,6 +96,13 @@ impl GLAreaImpl for Render {
             .borrow()
             .draw(canvas, self.scale.borrow().clone(), dpi);
 
+        self.foreground.borrow().draw(
+            canvas,
+            self.scale.borrow().clone(),
+            dpi,
+            self.mapper.borrow(),
+        );
+
         canvas.flush();
 
         true

src/render/mod.rs

@@ -2,13 +2,16 @@ mod background;
 mod foreground;
 
 mod imp;
+use std::fmt::Debug;
+
 use crate::data::{mapper::Mapper, RadarData2d};
 
 pub use self::background::{BackgroundConfig, BackgroundWidget};
 pub use self::foreground::{ForegroundConfig, ForegroundWidget};
+use crate::data::DownSampleMeth;
 pub use glib::subclass::prelude::*;
-use ndarray::{self, s, Axis, Dimension, Ix2, Zip};
+use ndarray::{self, s, Array2, Axis, Dimension, Ix2, Zip};
-use num_traits::Num;
+use num_traits::{AsPrimitive, FromPrimitive, Num};
 use proj::ProjError;
 
 pub(super) type WindowCoord = (f32, f32);
@@ -39,17 +42,30 @@ impl Render {
 
     pub(super) fn load_data_2d<T, Raw>(&self, data: RadarData2d<T, Raw>) -> Result<(), ProjError>
     where
-        T: Num + Clone + PartialEq + PartialOrd,
+        T: Num
+            + Clone
+            + PartialEq
+            + PartialOrd
+            + AsPrimitive<i8>
+            + AsPrimitive<f64>
+            + Debug
+            + FromPrimitive,
         Raw: ndarray::Data<Elem = T> + Clone + ndarray::RawDataClone,
     {
         assert!(data.dim1.shape().len() == data.dim2.shape().len());
 
+        data.downsample((801 / 2, 947 / 2), DownSampleMeth::VAR);
+
         let meshed = data.map_by_fn(|(x, y)| Ok((x, y)))?;
+        let d = data.data.to_owned().map(|x| x.as_());
 
         self.imp()
             .foreground
             .borrow_mut()
             .set_dims((meshed.dim1, meshed.dim2));
+
+        self.imp().foreground.borrow_mut().set_data(d);
+
         Ok(())
     }
 }

src/tree.rs

@@ -1,7 +1,7 @@
-use gtk::gdk::RGBA;
+use femtovg::Color;
 use num_traits::Num;
 
-pub fn get<T>(levels: &Vec<T>, colors: &Vec<RGBA>, v: T) -> RGBA
+pub fn get<T, V: Copy>(levels: &Vec<T>, colors: &Vec<V>, v: T) -> V
 where
     T: Num + PartialOrd + Copy,
 {