reliefImage.m

function outFiles = reliefDiff3D(noMocoInput, mocoInput, sliceFrac, mode, elevation, outDir)
%RELIEFDIFF3D  Relief of (moco - noMoco) across all 3 planes x 4 light angles.
%   For each anatomical plane (axial, coronal, sagittal) it takes a central
%   slice and renders the relief of the difference at four light azimuths
%   (0, 45, 90, 135 deg), so directional artefacts can't hide in a blind
%   orientation and you can inspect x, y and z. Saves one PNG per plane into
%   the moco file's sliding_L*_S* folder.
%
%   outFiles = reliefDiff3D(noMocoFile, mocoFile)
%   outFiles = reliefDiff3D(noMocoFile, mocoFile, sliceFrac, mode, elevation, outDir)
%   (arrays accepted, but pass PATHS so it can auto-route to the config folder)
%
%   INPUTS
%     noMocoInput, mocoInput : .mat paths (var auto-detected) or arrays.
%     sliceFrac : where to slice each axis, 0..1 (default 0.5 = centre).
%     mode      : 'directional' (default) or 'normal'.
%     elevation : light elevation for 'normal' (default pi/4).
%     outDir    : save folder (default = moco file's folder; pwd if array).
%
%   OUTPUT
%     outFiles : cell array of the 3 saved PNG paths (one per plane).
%
%   PLANES (for a [X Y Z] volume):
%     axial    = fixed Z (slice through dim 3) -> the view you've been using
%     coronal  = fixed Y (slice through dim 2)
%     sagittal = fixed X (slice through dim 1)

    if nargin < 3 || isempty(sliceFrac), sliceFrac = 0.5; end
    if nargin < 4 || isempty(mode),      mode = 'directional'; end
    if nargin < 5 || isempty(elevation), elevation = pi/4; end
    if nargin < 6, outDir = ''; end

    [A,~]        = loadVol(noMocoInput);
    [B,mocoPath] = loadVol(mocoInput);
    assert(isequal(size(A),size(B)),'Images must match in size (%s vs %s).', ...
        mat2str(size(A)), mat2str(size(B)));
    A = abs(double(A)); B = abs(double(B));
    A = A/prctile(A(:),99); B = B/prctile(B(:),99);   % common scale
    D = B - A;                                         % the difference volume
    sz = size(D);

    if isempty(outDir)
        if ~isempty(mocoPath), outDir = fileparts(mocoPath); else, outDir = pwd; end
    end
    if ~exist(outDir,'dir'), mkdir(outDir); end
    LSstr = parseLS(mocoPath); if isempty(LSstr), LSstr='config'; end

    angles = [0 pi/4 pi/2 3*pi/4];
    angLabels = {'0','45','90','135'};

    % plane definitions: name, slice index, how to extract the 2-D slice
    planes(1).name='axial';    planes(1).idx=round(sliceFrac*sz(3));
    planes(1).get=@(V,i) squeeze(V(:,:,i));
    planes(2).name='coronal';  planes(2).idx=round(sliceFrac*sz(2));
    planes(2).get=@(V,i) squeeze(V(:,i,:));
    planes(3).name='sagittal'; planes(3).idx=round(sliceFrac*sz(1));
    planes(3).get=@(V,i) squeeze(V(i,:,:));

    outFiles = {};
    for pp = 1:numel(planes)
        P = planes(pp);
        Dslice = P.get(D, clampIdx(P.idx, planeDim(pp), sz));

        fig = figure('Visible','off','Color','w','Position',[40 40 1300 360]);
        tl = tiledlayout(fig,1,4,'Padding','compact','TileSpacing','compact');
        for a = 1:numel(angles)
            r = relief2D(Dslice, angles(a), mode, elevation);
            nexttile(tl);
            if strcmpi(mode,'directional')
                lim = prctile(abs(r(:)),99)+eps;
                imagesc(r,[-lim lim]);
            else
                imagesc(r,[0 1]);
            end
            axis image off; colormap gray;
            title(sprintf('light %s\\circ', angLabels{a}));
        end
        sgtitle(sprintf('Relief of (moco - noMoco) — %s plane — %s', P.name, LSstr), ...
            'Interpreter','tex');

        f = fullfile(outDir, sprintf('reliefDiff_%s_%s_allAngles.png', P.name, LSstr));
        exportgraphics(fig, f, 'Resolution', 200);
        close(fig);
        outFiles{end+1} = f; %#ok<AGROW>
        fprintf('Saved %s-plane relief:\n%s\n', P.name, f);
    end
end

% ================= helpers =================
function d = planeDim(pp)
    % which volume dimension is held fixed for each plane
    switch pp, case 1, d=3; case 2, d=2; case 3, d=1; end
end
function i = clampIdx(i, dim, sz)
    i = max(1, min(i, sz(dim)));
end

function [V,pth] = loadVol(f)
    pth='';
    if isnumeric(f), V=f; if iscell(V), V=V{1}; end; return; end
    pth=f; S=load(f); fn=fieldnames(S); V=[];
    prefer={'x_moco','x_nomoco','x0','x_refined_moco','x'};
    for k=1:numel(prefer), if isfield(S,prefer{k}), V=S.(prefer{k}); break; end; end
    if isempty(V)
        for i=1:numel(fn)
            v=S.(fn{i}); if iscell(v)&&~isempty(v), v=v{1}; end
            if isnumeric(v)&&ndims(v)==3&&min(size(v))>8, V=v; break; end
        end
    end
    if iscell(V), V=V{1}; end
    assert(~isempty(V),'No 3-D array found in %s.',f);
end

function ls = parseLS(pth)
    ls='';
    if isempty(pth), return; end
    tok=regexp(pth,'(L\d+(?:\.\d+)?_S\d+(?:\.\d+)?)','tokens','once');
    if ~isempty(tok), ls=tok{1}; end
end

function relief = relief2D(img, theta, mode, elevation)
    x=double(img);                    % signed difference -> no abs
    gx=Dx1(x); gy=Dy1(x);
    switch lower(mode)
        case 'directional'
            relief=cos(theta)*gx+sin(theta)*gy;
        case 'normal'
            nrm=sqrt(gx.^2+gy.^2+1);
            nX=-gx./nrm; nY=-gy./nrm; nZ=1./nrm;
            lx=cos(elevation)*cos(theta); ly=cos(elevation)*sin(theta); lz=sin(elevation);
            relief=max(nX*lx+nY*ly+nZ*lz,0);
        otherwise, error('mode must be directional or normal.');
    end
end

function f = Dx1(Mat), [m,n]=size(Mat); f=(Mat([2:m m],1:n)-Mat([1 1:m-1],1:n))./2; end
function f = Dy1(Mat), [m,n]=size(Mat); f=(Mat(1:m,[2:n n])-Mat(1:m,[1 1:n-1]))./2; end

x_moco = 'I:\yannic.delisle\scripts\sub-10\sliding_L16.0_S3.2\x_moco_mathilda_smoothedMotion_from_x_cs_tres3.2s_Nu120x120x120_delta2_sliding.mat'; 
noMocoFile = 'I:\yannic.delisle\scripts\sub-10\baseline_nomoco\x_nomoco_1bin_mathilda_manual_finalNu240.mat';
T = load(mocoFile, 'x_moco');
moco = T.x_moco;
S = load(noMocoFile);              
x = S.x_nomoco;

slice = abs(moco(:,:,120));

reliefDiff3D(noMocoFile, x_moco);              % 45° light, saved to the L5.5_S5.5 folder