% Function part of NanoLocz GUI and NanoLocz-lib (2025).
%
% open_h5jpk - Open JPK AFM image data stored in HDF5 (.h5-jpk) format.
%
%   This function reads an HDF5 file generated by JPK AFM software, extracts
%   the available imaging channels, associated scan parameters, and the
%   requested image data, and returns the image stack and metadata.
%
%   Syntax:
%       [im, s] = open_h5jpk(f, ch)
%
%   Inputs:
%       f  - Full file path or filename of the HDF5 (.h5-jpk) file.
%       ch - Desired channel name as a string (e.g., 'Height-retrace', 'Height-trace').
%
%   Outputs:
%       im - 3D image array of size (xPixel, yPixel, number of frames). Pixel values
%            are scaled according to the selected channel and converted to nanometers
%            for Height channels.
%
%       s  - Structure containing scan metadata:
%            s.channels       - Cell array of available channels (with '-trace' or '-retrace' labels)
%            s.channel        - The actual channel used for image extraction
%            s.LineRate       - Scan line rate (Hz)
%            s.xPixel, s.yPixel - Number of pixels in X and Y directions
%            s.xsize, s.ysize - Scan size in X and Y (meters)
%            s.numberofFrames - Number of image frames
%
%   Notes:
%   - If the requested channel 'ch' is not found, the function defaults to 'Height-trace' 
%   - Pixel values are scaled using the stored multiplier and offset. For height channels, 
%     the result is converted to nanometers.
%
%   Example:
%       [im, s] = open_h5jpk('scan.jpk-force-map', 'Height-retrace');
%       imagesc(im(:,:,1)); axis image; colorbar;
%
function [im, s ] = open_h5jpk(f, ch)
info = h5info(f);
n_ch = numel(info.Groups.Groups)-1;
for i = 1:n_ch
    attrs = info.Groups.Groups(i).Attributes;
    names = {attrs.Name};
    fancy_idx = strcmp(names, 'channel.fancy-name');
    fancyName = attrs(fancy_idx).Value;
    type_idx = strcmp(names, 'retrace');
    channelType = attrs(type_idx).Value;

    if strcmp(channelType, 'true')
        s.channels{i} = [fancyName, '-retrace'];
    else
        s.channels{i} = [fancyName, '-trace'];
    end
end
attrs = info.Groups.Attributes;
names = {attrs.Name};

idx = strcmp(names, 'timing-settings.scanRate');
s.LineRate = attrs(idx).Value;
idx = strcmp(names, 'position-pattern.grid.ilength');
s.xPixel = attrs(idx).Value;
idx = strcmp(names, 'position-pattern.grid.jlength');
s.yPixel = attrs(idx).Value;
idx = strcmp(names, 'position-pattern.grid.ulength');
s.xsize = attrs(idx).Value;
idx = strcmp(names, 'position-pattern.grid.vlength');
s.ysize = attrs(idx).Value;

found_ch = find(strcmp(s.channels, ch));

if  any(found_ch)
    s.channel = ch;
else
    s.channel = 'Height-trace';
    ch = 'Height-retrace';
    
end

id_ch = find(strcmp(s.channels, ch));
attrs = info.Groups.Groups(id_ch).Attributes;
names = {attrs.Name};

% Try net-encoder first
% Default to raw encoder scaling
offs = attrs(strcmp(names, 'encoder.scaling.offset')).Value;
multi = attrs(strcmp(names, 'encoder.scaling.multiplier')).Value;
net = 0;

% Try net-encoder first
if any(strcmp(names, 'net-encoder.scaling.multiplier'))
    offs = attrs(strcmp(names, 'net-encoder.scaling.offset')).Value;
    multi = attrs(strcmp(names, 'net-encoder.scaling.multiplier')).Value;
    net = 1;

% Try calibrated path if available
elseif any(strcmp(names, 'conversion.calibrated.scaling.multiplier'))
    offs = attrs(strcmp(names, 'conversion.calibrated.scaling.offset')).Value;
    multi = attrs(strcmp(names, 'conversion.calibrated.scaling.multiplier')).Value;

    multi_nominal = attrs(strcmp(names, 'conversion.nominal.scaling.multiplier')).Value;
    offs_nominal = attrs(strcmp(names, 'conversion.nominal.scaling.offset')).Value;
    multi_encoder = attrs(strcmp(names, 'encoder.scaling.multiplier')).Value;
    offs_encoder = attrs(strcmp(names, 'encoder.scaling.offset')).Value;

    multi = multi * multi_nominal * multi_encoder;
    offs = offs_encoder * multi_nominal * multi + offs_nominal * multi + offs;
    net = 1;

% Try amplitude-specific conversion (e.g. distanceamplitude)
elseif any(strcmp(names, 'conversion.distanceamplitude.scaling.multiplier'))
    offs = attrs(strcmp(names, 'conversion.distanceamplitude.scaling.offset')).Value;
    multi = attrs(strcmp(names, 'conversion.distanceamplitude.scaling.multiplier')).Value;

    multi_amp = attrs(strcmp(names, 'conversion.voltsamplitude.scaling.multiplier')).Value;
    offs_amp = attrs(strcmp(names, 'conversion.voltsamplitude.scaling.offset')).Value;

    multi_encoder = attrs(strcmp(names, 'encoder.scaling.multiplier')).Value;
    offs_encoder = attrs(strcmp(names, 'encoder.scaling.offset')).Value;

    multi = multi * multi_amp * multi_encoder;
    offs = offs_encoder * multi_amp * multi + offs_amp * multi + offs;
    net = 1;
end

image_group = info.Groups.Groups(id_ch);
dataset_name = image_group.Datasets(1).Name;
full_path = [image_group.Name '/' dataset_name];
img_data = h5read(f, full_path);

s.numberofFrames = size(img_data, 1)-1;
im = zeros(s.xPixel, s.yPixel, s.numberofFrames);

for i = 1:s.numberofFrames
    im(:, :, i) = reshape(img_data(i, :), [s.xPixel, s.yPixel])';
end

im = (flip(im)*multi)+offs;
switch s.channel
    case {'Height-retrace','Height-trace'}
        if net==1
            im = im*1e9;
        end
    otherwise
end
end