Implement pyraview.readFile in MATLAB and Python

docs/API.md

@@ -1,6 +1,8 @@
 # Pyraview API Reference
 
-## C API (`src/c/pyraview.c`, `include/pyraview_header.h`)
+## C/C++ API (`src/c/pyraview.cpp`, `include/pyraview_header.h`)
+
+The core implementation is written in C++11 for efficient multi-threading, but exposes a C-compatible API for easy integration with other languages.
 
 ### `pyraview_process_chunk`
 Processes a chunk of raw data and updates decimation pyramids.
@@ -26,15 +28,15 @@ Arguments:
 - `levelSteps`: Pointer to array of decimation factors (e.g., `[100, 10, 10]`).
 - `numLevels`: Number of elements in `levelSteps`.
 - `nativeRate`: Original sampling rate (Hz).
-- `numThreads`: Number of OpenMP threads (0 for auto).
+- `numThreads`: Number of worker threads (0 for auto).
 
 Returns:
 - 0 on success.
 - Negative values on error (e.g., -2 mismatch, -1 I/O error).
 
 ---
 
-## Python API (`src/python/pyraview.py`)
+## Python API (`src/python/pyraview/__init__.py`)
 
 ### `process_chunk(data, file_prefix, level_steps, native_rate, append=False, layout='SxC', num_threads=0)`
 Wrapper for the C function.
@@ -52,6 +54,19 @@ Arguments:
 Returns:
 - 0 on success. Raises `RuntimeError` on failure.
 
+### `read_file(filename, s0, s1)`
+Reads a specific range of samples from a level file.
+
+Arguments:
+- `filename`: String path to the file.
+- `s0`: Start sample index (int or float). Use `float('-inf')` for beginning.
+- `s1`: End sample index (int or float). Use `float('inf')` for end.
+
+Returns:
+- Numpy array of shape `(Samples, Channels, 2)`.
+    - `[:, :, 0]`: Minimum values.
+    - `[:, :, 1]`: Maximum values.
+
 ---
 
 ## Matlab API (`src/matlab/pyraview_mex.c`)
@@ -68,3 +83,16 @@ Arguments:
 
 Returns:
 - `status`: 0 on success. Throws error on failure.
+
+### `D = pyraview.readFile(filename, s0, s1)`
+Reads a specific range of samples from a level file.
+
+Arguments:
+- `filename`: String path to the `.bin` level file.
+- `s0`: Start sample index (0-based). Can be `-Inf`.
+- `s1`: End sample index (0-based). Can be `Inf`.
+
+Returns:
+- `D`: A 3D matrix of size `(Samples x Channels x 2)`.
+    - `D(:, :, 1)`: Minimum values.
+    - `D(:, :, 2)`: Maximum values.

src/matlab/+pyraview/readFile.m

@@ -0,0 +1,139 @@
+function d = readFile(filename, s0, s1)
+%READFILE Reads a specific range of samples from a Pyraview level file.
+%
+%   D = pyraview.readFile(FILENAME, S0, S1) reads data from the file specified
+%   by FILENAME, starting at sample index S0 and ending at sample index S1
+%   (inclusive, 0-based indexing).
+%
+%   Inputs:
+%       FILENAME - String or character vector specifying the path to the Pyraview level file.
+%       S0       - Scalar numeric. The starting sample index (0-based).
+%                  Can be -Inf to indicate the beginning of the file.
+%       S1       - Scalar numeric. The ending sample index (0-based).
+%                  Can be Inf to indicate the end of the file.
+%
+%   Outputs:
+%       D        - A 3-D matrix of size (Samples x Channels x 2).
+%                  The data type matches the file's internal data type.
+%                  D(:, :, 1) contains the minimum values for each sample.
+%                  D(:, :, 2) contains the maximum values for each sample.
+%
+%   Example:
+%       % Read samples 100 to 200 from 'data_L1.bin'
+%       d = pyraview.readFile('data_L1.bin', 100, 200);
+%
+%       % Read from the beginning to sample 500
+%       d = pyraview.readFile('data_L1.bin', -Inf, 500);
+%
+%       % Read from sample 1000 to the end
+%       d = pyraview.readFile('data_L1.bin', 1000, Inf);
+%
+%   Notes:
+%       - Indices are clamped to the file's valid sample range.
+%       - Returns an empty array if the requested range is invalid or empty.
+%       - Pyraview files use a planar layout where channels are stored contiguously.
+%       - Each "sample" in a level file consists of a Min/Max pair.
+
+    if ~isfile(filename)
+        error('Pyraview:FileNotFound', 'File not found: %s', filename);
+    end
+
+    % Read header
+    try
+        h = pyraview.pyraview_get_header_mex(filename);
+    catch e
+        error('Pyraview:HeaderError', 'Failed to read header: %s', e.message);
+    end
+
+    % Determine data type and item size
+    switch h.dataType
+        case 0, precision = 'int8'; itemSize = 1;
+        case 1, precision = 'uint8'; itemSize = 1;
+        case 2, precision = 'int16'; itemSize = 2;
+        case 3, precision = 'uint16'; itemSize = 2;
+        case 4, precision = 'int32'; itemSize = 4;
+        case 5, precision = 'uint32'; itemSize = 4;
+        case 6, precision = 'int64'; itemSize = 8;
+        case 7, precision = 'uint64'; itemSize = 8;
+        case 8, precision = 'single'; itemSize = 4;
+        case 9, precision = 'double'; itemSize = 8;
+        otherwise
+            error('Pyraview:UnknownType', 'Unknown data type code: %d', h.dataType);
+    end
+
+    % Determine total samples
+    dDir = dir(filename);
+    fileSize = dDir.bytes;
+    headerSize = 1024;
+    dataArea = fileSize - headerSize;
+
+    numChannels = double(h.channelCount);
+    frameSize = numChannels * 2 * itemSize;
+    totalSamples = floor(dataArea / frameSize);
+
+    % Handle s0 and s1
+    if isinf(s0) && s0 < 0
+        startSample = 0;
+    else
+        startSample = s0;
+    end
+
+    if isinf(s1) && s1 > 0
+        endSample = totalSamples - 1;
+    else
+        endSample = s1;
+    end
+
+    % Validate indices
+    if startSample < 0
+        startSample = 0;
+    end
+
+    if endSample >= totalSamples
+        endSample = totalSamples - 1;
+    end
+
+    if startSample > endSample
+        d = zeros(0, numChannels, 2, precision);
+        return;
+    end
+
+    numSamplesToRead = endSample - startSample + 1;
+
+    % Allocate output
+    d = zeros(numSamplesToRead, numChannels, 2, precision);
+
+    % Open file
+    f = fopen(filename, 'rb');
+    if f == -1
+        error('Pyraview:FileOpenError', 'Could not open file: %s', filename);
+    end
+
+    cleanupObj = onCleanup(@() fclose(f));
+
+    for ch = 1:numChannels
+        % Calculate offset
+        % Planar layout: [Header] [Ch1 Data] [Ch2 Data] ...
+        % Ch Data size = totalSamples * 2 * itemSize
+
+        chStartOffset = headerSize + (ch-1) * totalSamples * 2 * itemSize;
+        readOffset = chStartOffset + startSample * 2 * itemSize;
+
+        fseek(f, readOffset, 'bof');
+
+        % Read min/max pairs
+        % precision needs to be char for fread
+        raw = fread(f, numSamplesToRead * 2, ['*' char(precision)]);
+
+        if length(raw) < numSamplesToRead * 2
+            warning('Pyraview:ShortRead', 'Short read on channel %d. Expected %d, got %d.', ch, numSamplesToRead*2, length(raw));
+            % Fill what we got
+            nRead = floor(length(raw)/2);
+            d(1:nRead, ch, 1) = raw(1:2:2*nRead);
+            d(1:nRead, ch, 2) = raw(2:2:2*nRead);
+        else
+            d(:, ch, 1) = raw(1:2:end);
+            d(:, ch, 2) = raw(2:2:end);
+        end
+    end
+end

src/matlab/tests/+pyraview/+unittest/test_readfile.m

@@ -0,0 +1,127 @@
+classdef test_readfile < matlab.unittest.TestCase
+
+    properties
+        TempFilename
+        CleanupObj
+        NumChannels = 2
+        NumSamples = 10
+        DataType = 'int16'
+        ItemSize = 2
+        DataCh1
+        DataCh2
+    end
+
+    methods(TestClassSetup)
+        function setupPaths(testCase)
+            % Verify MEX file exists
+            [~, mexName] = fileparts('pyraview');
+            mexExt = mexext;
+            fullMexPath = fullfile(pwd, 'src', 'matlab', ['pyraview.' mexExt]);
+
+            if ~exist(fullMexPath, 'file')
+                % Try relative to where this file is
+                currentFileDir = fileparts(mfilename('fullpath'));
+                fullMexPath = fullfile(currentFileDir, ['pyraview.' mexExt]);
+                if ~exist(fullMexPath, 'file')
+                     warning('MEX file not found at expected location.');
+                end
+                addpath(currentFileDir);
+                addpath(fullfile(currentFileDir, '+pyraview'));
+            else
+                addpath(fileparts(fullMexPath));
+                addpath(fullfile(fileparts(fullMexPath), '+pyraview'));
+            end
+        end
+    end
+
+    methods(TestMethodSetup)
+        function createTestFile(testCase)
+            testCase.TempFilename = [tempname '.bin'];
+            testCase.CleanupObj = onCleanup(@() delete(testCase.TempFilename));
+
+            % Generate data
+            testCase.DataCh1 = zeros(testCase.NumSamples * 2, 1, testCase.DataType);
+            for i = 1:testCase.NumSamples
+                testCase.DataCh1((i-1)*2 + 1) = (i-1)*2;     % Min
+                testCase.DataCh1((i-1)*2 + 2) = (i-1)*2 + 1; % Max
+            end
+
+            testCase.DataCh2 = zeros(testCase.NumSamples * 2, 1, testCase.DataType);
+            for i = 1:testCase.NumSamples
+                testCase.DataCh2((i-1)*2 + 1) = 100 + (i-1)*2;     % Min
+                testCase.DataCh2((i-1)*2 + 2) = 100 + (i-1)*2 + 1; % Max
+            end
+
+            % Write file
+            fid = fopen(testCase.TempFilename, 'wb');
+
+            % Header
+            fwrite(fid, 'PYRA', 'char');
+            fwrite(fid, 1, 'uint32');
+            fwrite(fid, 2, 'uint32'); % int16
+            fwrite(fid, testCase.NumChannels, 'uint32');
+            fwrite(fid, 1000.0, 'double');
+            fwrite(fid, 1000.0, 'double');
+            fwrite(fid, 0.0, 'double');
+            fwrite(fid, 1, 'uint32');
+            fwrite(fid, zeros(980, 1, 'uint8'), 'uint8');
+
+            % Data (Planar)
+            fwrite(fid, testCase.DataCh1, testCase.DataType);
+            fwrite(fid, testCase.DataCh2, testCase.DataType);
+
+            fclose(fid);
+        end
+    end
+
+    methods(TestMethodTeardown)
+        function deleteTestFile(testCase)
+            delete(testCase.CleanupObj);
+        end
+    end
+
+    methods(Test)
+        function testReadFull(testCase)
+            d = pyraview.readFile(testCase.TempFilename, 0, testCase.NumSamples-1);
+
+            testCase.verifyEqual(size(d), [testCase.NumSamples, testCase.NumChannels, 2]);
+            testCase.verifyEqual(d(:, 1, 1), testCase.DataCh1(1:2:end)); % Ch1 Mins
+            testCase.verifyEqual(d(:, 1, 2), testCase.DataCh1(2:2:end)); % Ch1 Maxs
+            testCase.verifyEqual(d(:, 2, 1), testCase.DataCh2(1:2:end)); % Ch2 Mins
+            testCase.verifyEqual(d(:, 2, 2), testCase.DataCh2(2:2:end)); % Ch2 Maxs
+        end
+
+        function testReadPartial(testCase)
+            s0 = 2; s1 = 4;
+            d = pyraview.readFile(testCase.TempFilename, s0, s1);
+
+            numRead = s1 - s0 + 1;
+            testCase.verifyEqual(size(d), [numRead, testCase.NumChannels, 2]);
+
+            % Check Ch1 Mins
+            expectedMins = testCase.DataCh1(1:2:end);
+            expectedMins = expectedMins(s0+1 : s1+1); % Matlab 1-based indexing
+            testCase.verifyEqual(d(:, 1, 1), expectedMins);
+        end
+
+        function testReadInf(testCase)
+            d = pyraview.readFile(testCase.TempFilename, -Inf, Inf);
+
+            testCase.verifyEqual(size(d), [testCase.NumSamples, testCase.NumChannels, 2]);
+            testCase.verifyEqual(d(:, 1, 1), testCase.DataCh1(1:2:end));
+        end
+
+        function testReadEmpty(testCase)
+            d = pyraview.readFile(testCase.TempFilename, 5, 4);
+            testCase.verifyEqual(size(d), [0, testCase.NumChannels, 2]);
+        end
+
+        function testReadOutOfBounds(testCase)
+            % Request beyond end, should be clamped
+            d = pyraview.readFile(testCase.TempFilename, testCase.NumSamples-2, testCase.NumSamples+100);
+
+            numRead = 2; % Samples 8 and 9 (0-based)
+            testCase.verifyEqual(size(d, 1), numRead);
+        end
+    end
+end