function overlay_grid(nii_file,structural_dir)
    %% 
    % 18-4-2025 Paul, MRSHub
    % Only for research use.
    % eg nii_file = 'G:\spectro_project\csi\csi.nii'
    % eg structural_dir = 'G:\spectro_project\im\'  (contains *.nii.gz')
    % Make sure pyenv points to python.exe

    %% 1. Load CSI
    pyenv('Version','C:\...\Local\Programs\Python\Python310\python.exe');
    nib     = py.importlib.import_module('nibabel');

    img    = nib.load(nii_file);
    h      = img.header;
    shp    = cellfun(@double, cell(h.get_data_shape()));  % [nx,ny,nz]
    A      = double(img.affine);

    % decide left/right flip from CSI affine x‑axis direction
    xd       = A(1:3,1) / norm(A(1:3,1));   % direction cosine of i‑axis
    doLRflip = (xd(1) < 0);                 % true if i→ decreases world X

    [X,Y]  = ndgrid(0:shp(1)-1, 0:shp(2)-1);
    Vc     = [X(:), Y(:), zeros(numel(X),1), ones(numel(X),1)]';
    worldC = A * Vc;     worldC = worldC(1:3,:);

    U        = img.dataobj.get_unscaled();
    Rdat     = double(py.numpy.array(U.real));
    Idat     = double(py.numpy.array(U.imag));
    csi_data = reshape(complex(Rdat,Idat), shp);

%%

 % Extract voxel dimensions.
 header=h;
    try
        voxel_dims = header.get_zooms();
        voxel_dims_matlab = cellfun(@double, cell(voxel_dims));
    catch
        voxel_dims_matlab = [];
    end
    
    % Extract header extension information.
    extList = header.extensions;
    if ~isempty(extList)
        extListCell = cell(extList);
        extObj = extListCell{1};
        extContent = char(extObj.get_content());
        if startsWith(extContent, "b'")
            extContent = extContent(3:end);
        end
        if endsWith(extContent, "'")
            extContent = extContent(1:end-1);
        end
        headerExt = jsondecode(extContent);
    else
        headerExt = [];
    end

    % Create header structure.
    spectroHeader = struct();
    spectroHeader.voxel_dims = voxel_dims_matlab;
    spectroHeader.header_extension = headerExt;
    assignin('base', 'spectroHeader', spectroHeader);

    %% 2. Load structural & invert affine
    Sfiles = dir(fullfile(structural_dir,'*.nii.gz'));
    volF   = fullfile(Sfiles(1).folder, Sfiles(1).name);
    Sinfo  = niftiinfo(volF);
    Svol   = niftiread(volF);
    T      = Sinfo.Transform.T';    % world ← voxel
    invT   = inv(T);                % voxel ← world
    pix    = Sinfo.PixelDimensions; % [dx,dy,dz]
    zd     = cellfun(@double, cell(h.get_zooms()));  % CSI [dx,dy,dz]

    %% 2.1 Precompute CSI centres & swap X↔Y
    tmp      = invT * A * Vc;
    newCoord = tmp(1:3,:);            
    newCoord = newCoord([2 1 3],:);   % (1)=cols, (2)=rows

    %% 3. CSI→pixel sizing (for rotated display)
    voxelW = zd(1)/pix(1);   % original X‑size in pixels
    voxelH = zd(2)/pix(2);   % original Y‑size in pixels
    % after +90° rotation, width→orig Y, height→orig X
    boxW  = voxelH;
    boxH  = voxelW;
    halfW = boxW/2;
    halfH = boxH/2;

    %% 4. GUI
    [nR,nC,nS] = size(Svol);
    f = figure('Name','CSI Overlay');
    ax = axes(f);

    if nS>1
      st = 1/(nS-1);
      sliderStep = [st, min(10/(nS-1),1)];
    else
      sliderStep = [1,1];
    end
    uicontrol(f,'Style','slider','Units','norm', ...
      'Pos',[.25 .02 .5 .05],'Min',1,'Max',nS,'Value',1, ...
      'SliderStep',sliderStep,'Callback',@(s,~) updateSlice(round(s.Value)));
    txt = uicontrol(f,'Style','text','Units','norm','Pos',[.76 .02 .1 .05], ...
      'String','Slice: 1');
    uicontrol(f,'Style','togglebutton','Units','norm','Pos',[.3 .9 .15 .05], ...
      'String','Voxel Select Mode','Callback',@toggleSelect);
    uicontrol(f,'Style','pushbutton','Units','norm','Pos',[.1 .9 .15 .05], ...
      'String','Extract Spectrum','Callback',@extractSpectrum);
    uicontrol(f,'Style','pushbutton','Units','norm','Pos',[.5 .9 .15 .05], ...
      'String','Average Selected','Callback',@averageSelected);

    pixelCenters=[]; pixelValidIdx=[]; selected=[]; currentSlice=1;
    updateSlice(1);

    %% updateSlice
    function updateSlice(k)
      currentSlice = k;
      I2 = double(Svol(:,:,k));
      if max(I2(:))>0, I2 = I2/max(I2(:)); end

      % rotate structural by +90° CCW
      R1 = rot90(I2,1);
      % optionally mirror L↔R
      if doLRflip
        R1 = fliplr(R1);
      end
      imshow(R1,'Parent',ax); hold(ax,'on');
      txt.String = ['Slice: ' num2str(k)];

      % slice‑mask in world space
      base = T*[0;0;k-1;1];
      p2   = T*[1;0;k-1;1];
      p3   = T*[0;1;k-1;1];
      nrm  = cross(p2(1:3)-base(1:3), p3(1:3)-base(1:3));
      nrm  = nrm / norm(nrm);
      d    = worldC - base(1:3);
      mask = abs(nrm'*d) < (zd(3)/2);
      idx  = find(mask);

      % map CSI centres through rotation (no pre‑swap here)
      coords = newCoord(:,mask);
      c0     = coords(1,:);  % original col
      r0     = coords(2,:);  % original row

      % after +90°: newCol = r0, newRow = nC - c0 +1
      newC = r0;
      newR = nC - c0 + 1;
      % apply L↔R flip if needed
      if doLRflip
        newC = size(R1,2) - newC + 1;
      end

      % convert to 1‑based pixel centers + snap to integers
      newC = newC - 1;
      newR = newR - 1;

      pixelCenters  = [newC; newR];
      pixelValidIdx = idx;

      % draw each box with swapped dims
      for i=1:numel(newC)
        sel = ismember(idx(i),selected);
        clr = sel*[1 1 0] + ~sel*[0 1 0];
        lw  = sel*2 + 1;

        rectangle('Position', ...
          [ newC(i)-halfH, newR(i)-halfW, boxH, boxW ], ...
          'EdgeColor',clr,'LineWidth',lw);

      end
      hold(ax,'off');
    end

    %% voxel‑select toggle & pickVoxel
    function toggleSelect(src,~)
      if src.Value
        f.WindowButtonDownFcn = @pickVoxel;
        src.String = 'Selection ON';
      else
        f.WindowButtonDownFcn = [];
        src.String = 'Voxel Select Mode';
      end
    end
    function pickVoxel(~,~)
      cp = ax.CurrentPoint(1,1:2);
      d  = sqrt(sum((pixelCenters'-cp).^2,2));
      [md,i] = min(d);
      if md < max(halfW,halfH)*1.5
        li = pixelValidIdx(i);
        if ismember(li,selected)
          selected(selected==li) = [];
        else
          selected(end+1) = li;
        end
        updateSlice(currentSlice);
      end
    end

    %% extractSpectrum
    function extractSpectrum(~,~)
      disp('Click on a CSI voxel…');
      [x,y] = ginput(1);
      updateSlice(currentSlice);
      d  = sqrt(sum((pixelCenters'-[x y]).^2,2));
      [~,i] = min(d);
      idx0   = pixelValidIdx(i);
      [rr,cc]= ind2sub([shp(1),shp(2)],idx0);
      rawSpec= squeeze(csi_data(rr,cc,1,:,1));

      assignin('base', 'extractedSpectrum', rawSpec);

      figure;
      subplot(2,1,1);
      plot(real(rawSpec)); hold on; plot(imag(rawSpec)); hold off;
      xlabel('Point'); ylabel('Signal');
      title(sprintf('Raw Spectrum (%d,%d)',rr,cc));
      legend('Real','Imag');
      subplot(2,1,2);
      plot(abs(fftshift(fft(rawSpec))),'LineWidth',1.5);
      xlabel('Freq Index'); ylabel('Magnitude');
      title('FFT Magnitude');
    end

    %% averageSelected
    function averageSelected(~,~)
      if isempty(selected), disp('No voxels selected.'), return, end
      N     = numel(selected);
      specs = cell(N,1); minL = inf;
      for i=1:N
        [rr,cc] = ind2sub([shp(1),shp(2)],selected(i));
        s       = squeeze(csi_data(rr,cc,1,:,1));
        specs{i}= s;
        minL    = min(minL,numel(s));
      end
      M       = zeros(minL,N);
      for i=1:N, M(:,i)= specs{i}(1:minL); end
      avgSpec = mean(M,2);

      figure;
      subplot(2,1,1);
      plot(real(avgSpec)); hold on; plot(imag(avgSpec)); hold off;
      xlabel('Point'); ylabel('Signal');
      title('Averaged Time‑Domain Spectrum');
      legend('Real','Imag');
      subplot(2,1,2);
      plot(abs(fftshift(fft(avgSpec))),'LineWidth',1.5);
      xlabel('Freq Index'); ylabel('Magnitude');
      title('Averaged FFT Magnitude');

      assignin('base', 'averagedSpectrum', avgSpec);


    end
end
