wip arbitrary arrays

Author: normanrz

matlab/rust/wkw_load/src/lib.rs

@@ -4,7 +4,6 @@ extern crate zarrs;
 #[macro_use]
 extern crate wkw_mex;
 use wkw_mex::*;
-use wkwrap::{Box3, Mat, Vec3, VoxelType};
 use zarrs::array::data_type::DataType;
 use zarrs::array::Array;
 use zarrs::array_subset::ArraySubset;
@@ -40,55 +39,37 @@ mex_function!(nlhs, lhs, nrhs, rhs, {
     );
     let array = zarrs_result_to_str_error(Array::open(store.clone(), "/"))?;
 
-    let num_channels = array.shape()[0] as usize;
-    let is_multi_channel = num_channels > 1;
+    let array_shape = array.shape();
+    let ndim = array_shape.len();
+    let data_type = array.data_type();
+    let type_size = if let Some(type_size) = data_type.fixed_size() {
+        type_size
+    } else {
+        return Err("Unsupported data type".to_string());
+    };
 
     // build shape
-    let bbox = mx_array_to_wkwrap_box(rhs[1])?;
+    let (bbox_start, bbox_shape) = mx_array_to_bbox(rhs[1], ndim)?;
+
+    if bbox_start
+        .iter()
+        .zip(bbox_shape.iter())
+        .zip(array_shape.iter())
+        .any(|((bbox_min_x, bbox_shape_x), shape_x)| (*bbox_min_x + *bbox_shape_x) > *shape_x)
+    {
+        return Err(format!(
+            "Bounding box start={:?}, shape={:?} is out of bounds for array of shape={:?}.",
+            bbox_start, bbox_shape, array_shape
+        ));
+    }
+    println!("{:?} {:?}", bbox_start, bbox_shape);
     let subset = zarrs_result_to_str_error(ArraySubset::new_with_start_shape(
-        vec![
-            0,
-            bbox.min().x as u64,
-            bbox.min().y as u64,
-            bbox.min().z as u64,
-        ],
-        vec![
-            1,
-            bbox.width().x as u64,
-            bbox.width().y as u64,
-            bbox.width().z as u64,
-        ],
+        bbox_start.clone(),
+        bbox_shape.clone(),
     ))?;
 
-    let shape_arr = [
-        num_channels,
-        bbox.width().x as usize,
-        bbox.width().y as usize,
-        bbox.width().z as usize,
-    ];
-    let shape_slice = if is_multi_channel {
-        &shape_arr[0..]
-    } else {
-        &shape_arr[1..]
-    };
-
     // prepare allocation
-    let voxel_type = match array.data_type() {
-        DataType::UInt8 => VoxelType::U8,
-        DataType::UInt16 => VoxelType::U16,
-        DataType::UInt32 => VoxelType::U32,
-        DataType::UInt64 => VoxelType::U64,
-        DataType::Float32 => VoxelType::F32,
-        DataType::Float64 => VoxelType::F64,
-        DataType::Int8 => VoxelType::I8,
-        DataType::Int16 => VoxelType::I16,
-        DataType::Int32 => VoxelType::I32,
-        DataType::Int64 => VoxelType::I64,
-        _ => {
-            return Err("Unsupported data type".to_string());
-        }
-    };
-    let class = match array.data_type() {
+    let mat_class = match array.data_type() {
         DataType::UInt8 => MxClassId::Uint8,
         DataType::UInt16 => MxClassId::Uint16,
         DataType::UInt32 => MxClassId::Uint32,
@@ -106,24 +87,16 @@ mex_function!(nlhs, lhs, nrhs, rhs, {
 
     // read data
     let data_all = zarrs_result_to_str_error(array.retrieve_array_subset(&subset))?;
-    let mut buf = zarrs_result_to_str_error(data_all.into_fixed())?.into_owned(); // in c-order
-    let src_mat = Mat::new(
-        &mut buf,
-        bbox.width(),
-        voxel_type.size() * num_channels,
-        voxel_type,
-        true,
-    )?;
-    let arr = create_numeric_array(shape_slice, class, MxComplexity::Real)?;
-    let mut mat = mx_array_mut_to_wkwrap_mat(is_multi_channel, arr)?;
-
-    src_mat.copy_as_fortran_order(
-        &mut mat,
-        Box3::new(Vec3 { x: 0, y: 0, z: 0 }, bbox.width())?,
-    )?;
+    let zarr_buf = zarrs_result_to_str_error(data_all.into_fixed())?.into_owned(); // in c-order
+
+    println!("zarr_buf: {:?}", zarr_buf);
+
+    let mat_arr = create_numeric_array(&bbox_shape, mat_class, MxComplexity::Real)?;
+
+    copy_as_fortran_order(&zarr_buf, mat_arr, &bbox_shape, type_size)?;
 
     // set output
-    lhs[0] = arr;
+    lhs[0] = mat_arr;
 
     Ok(())
 });

matlab/rust/wkw_mex/src/util.rs

@@ -1,19 +1,19 @@
-use ::ffi::*;
+use ffi::*;
 
 use std;
-use std::slice;
 use std::ffi::CStr;
+use std::slice;
 
 pub type Result<T> = std::result::Result<T, String>;
 
 pub fn as_nat(f: f64) -> Result<u64> {
     if f <= 0.0 {
-        return Err("Input must be positive".to_string())
+        return Err("Input must be positive".to_string());
     }
 
     match f % 1.0 == 0.0 {
         true => Ok(f as u64),
-        false => Err("Input must be an integer".to_string())
+        false => Err("Input must be an integer".to_string()),
     }
 }
 
@@ -22,53 +22,57 @@ pub fn as_log2(f: f64) -> Result<u8> {
 
     match i & (i - 1) == 0 {
         true => Ok(i.trailing_zeros() as u8),
-        false => Err("Input must be a power of two".to_string())
+        false => Err("Input must be a power of two".to_string()),
     }
 }
 
 pub fn str_slice_to_mx_class_id(class_id: &str) -> Result<MxClassId> {
     match class_id {
-        "uint8"  => Ok(MxClassId::Uint8),
+        "uint8" => Ok(MxClassId::Uint8),
         "uint16" => Ok(MxClassId::Uint16),
         "uint32" => Ok(MxClassId::Uint32),
         "uint64" => Ok(MxClassId::Uint64),
         "single" => Ok(MxClassId::Single),
         "double" => Ok(MxClassId::Double),
-        "int8"   => Ok(MxClassId::Int8),
-        "int16"  => Ok(MxClassId::Int16),
-        "int32"  => Ok(MxClassId::Int32),
-        "int64"  => Ok(MxClassId::Int64),
-        _        => Err("Unknown MxClassId name".to_string())
+        "int8" => Ok(MxClassId::Int8),
+        "int16" => Ok(MxClassId::Int16),
+        "int32" => Ok(MxClassId::Int32),
+        "int64" => Ok(MxClassId::Int64),
+        _ => Err("Unknown MxClassId name".to_string()),
     }
 }
 
 pub fn mx_array_to_str<'a>(pm: MxArray) -> Result<&'a str> {
     let pm_ptr = unsafe { mxArrayToUTF8String(pm) };
 
     if pm_ptr.is_null() {
-        return Err("mxArrayToUTF8String returned null".to_string())
+        return Err("mxArrayToUTF8String returned null".to_string());
     }
 
     let pm_cstr = unsafe { CStr::from_ptr(pm_ptr) };
 
     match pm_cstr.to_str() {
         Ok(pm_str) => Ok(pm_str),
-        Err(_) => Err("mxArray contains invalid UTF-8 data".to_string())
+        Err(_) => Err("mxArray contains invalid UTF-8 data".to_string()),
     }
 }
 
 pub fn mx_array_to_f64_slice<'a>(pm: MxArray) -> Result<&'a [f64]> {
     unsafe {
-        if !mxIsDouble(pm) { return Err("MxArray is not of class \"double\"".to_string()) };
-        if mxIsComplex(pm) { return Err("MxArray is complex".to_string()) };
+        if !mxIsDouble(pm) {
+            return Err("MxArray is not of class \"double\"".to_string());
+        };
+        if mxIsComplex(pm) {
+            return Err("MxArray is complex".to_string());
+        };
     }
 
     let pm_numel = unsafe { mxGetNumberOfElements(pm) };
     let pm_ptr = unsafe { mxGetPr(pm) };
 
     match pm_ptr.is_null() {
         true => Err("MxArray does not contain real values".to_string()),
-        false => Ok(unsafe { slice::from_raw_parts(pm_ptr, pm_numel) })
+        false => Ok(unsafe { slice::from_raw_parts(pm_ptr, pm_numel) }),
     }
 }
 
@@ -77,7 +81,7 @@ pub fn mx_array_to_f64(pm: MxArray) -> Result<f64> {
 
     match pm_slice.len() {
         1 => Ok(pm_slice[0]),
-        _ => Err("MxArray contains an invalid number of doubles".to_string())
+        _ => Err("MxArray contains an invalid number of doubles".to_string()),
     }
 }
 
@@ -113,25 +117,26 @@ pub fn mx_array_size_to_usize_slice<'a>(pm: MxArray) -> &'a [usize] {
     let ndims = unsafe { mxGetNumberOfDimensions(pm) };
     let dims = unsafe { mxGetDimensions(pm) };
 
-    unsafe {
-        slice::from_raw_parts(dims, ndims as usize)
-    }
+    unsafe { slice::from_raw_parts(dims, ndims as usize) }
 }
 
 pub fn create_numeric_array(
-    dims: &[usize],
+    dims: &[u64],
     class: MxClassId,
-    complexity: MxComplexity
+    complexity: MxComplexity,
 ) -> Result<MxArrayMut> {
     let arr = unsafe {
         mxCreateNumericArray(
-            dims.len() as size_t, dims.as_ptr(),
-            class as c_int, complexity as c_int)
+            dims.len() as size_t,
+            dims.as_ptr() as *const usize,
+            class as c_int,
+            complexity as c_int,
+        )
     };
 
     match arr.is_null() {
         true => Err("Failed to create uninitialized numeric array".to_string()),
-        false => Ok(arr)
+        false => Ok(arr),
     }
 }
 
@@ -140,7 +145,7 @@ pub fn malloc(n: usize) -> Result<&'static mut [u8]> {
 
     match ptr.is_null() {
         true => Err("Failed to allocate memory".to_string()),
-        false => Ok(unsafe { slice::from_raw_parts_mut(ptr, n) })
+        false => Ok(unsafe { slice::from_raw_parts_mut(ptr, n) }),
     }
 }
 
@@ -156,3 +161,63 @@ pub fn die(msg: &str) {
     // die
     unsafe { mexErrMsgTxt(buf.as_ptr()) }
 }
+
+pub fn copy_as_fortran_order(
+    in_buf: &[u8],
+    out_arr: MxArrayMut,
+    shape: &[u64],
+    type_size: usize,
+) -> Result<()> {
+    let total_elems: usize = shape.iter().product::<u64>() as usize;
+    if in_buf.len() != total_elems * type_size {
+        return Err(format!(
+            "Length of input buffer does not match expected size {} != {}",
+            in_buf.len(),
+            total_elems,
+        ));
+    }
+
+    let result = mx_array_mut_to_u8_slice_mut(out_arr)?;
+    if result.len() != total_elems * type_size {
+        return Err(format!(
+            "Length of output array does not match expected size {} != {}",
+            result.len(),
+            total_elems,
+        ));
+    }
+
+    // Compute F-order (column-major) strides
+    let mut f_strides = vec![1u64; shape.len()];
+    for i in 1..shape.len() {
+        f_strides[i] = f_strides[i - 1] * shape[i - 1];
+    }
+
+    // Multi-dimensional index for C-order iteration
+    let mut idx = vec![0u64; shape.len()];
+
+    // Iterate over all elements in C-order (sequential read)
+    for elem_idx in 0..total_elems {
+        // Compute Fortran-order (column-major) offset
+        let f_offset_elems: u64 = idx.iter().zip(&f_strides).map(|(&i, &s)| i * s).sum();
+        let f_offset_bytes = f_offset_elems as usize * type_size;
+
+        // Sequential read from in_buf
+        let src_offset_bytes = elem_idx * type_size;
+        let src_slice = &in_buf[src_offset_bytes..src_offset_bytes + type_size];
+
+        // Scattered write to result
+        result[f_offset_bytes..f_offset_bytes + type_size].copy_from_slice(src_slice);
+
+        // Increment multi-dimensional index (C-order)
+        for d in (0..shape.len()).rev() {
+            idx[d] += 1;
+            if idx[d] < shape[d] {
+                break;
+            } else if d > 0 {
+                idx[d] = 0;
+            }
+        }
+    }
+
+    Ok(())
+}