% WRONG: Mixing coordinate systems
pixel = V.Voxels(100, 200, 50);  % Intrinsic indexing
world_point = [100, 200, 50];    % Treating indices as world coords!

% CORRECT: Explicit conversion (separate coordinate arrays)
V = medicalVolume('scan.nii');
i = 100; j = 200; k = 50;  % Voxel indices
[x, y, z] = intrinsicToWorld(V.VolumeGeometry, i, j, k);
% x, y, z are now in mm, in patient coordinate system
world_point = [x, y, z];

% Convert back (also requires separate arrays)
[row, col, slice] = worldToSubscript(V.VolumeGeometry, x, y, z);

% WRONG: Loses spatial referencing
data = niftiread('brain.nii');  % Just a 3D array, no spatial info

% CORRECT: Preserves geometry
V = medicalVolume('brain.nii');
% V.VolumeGeometry contains transformation matrix
% V.VoxelSpacing is in mm
% V.SpatialUnits documents the units

V = medicalVolume('scan.dcm');

% Check orientation
disp(V.Orientation);        % 'transverse', 'sagittal', or 'coronal'
disp(V.PlaneMapping);       % Maps dimensions to anatomical planes

% Extract slices in patient coordinates (not array indices!)
slice = extractSlice(V, 50, 'transverse');  % 50th transverse slice

% For axial browsing regardless of acquisition:
sliceViewer(V);  % Automatically handles orientation

% Get collection info for DICOM series
collection = dicomCollection('DICOM_folder');
disp(collection.Properties.VariableNames);

% Check modality, series description
info = dicominfo('slice001.dcm');
fprintf('Modality: %s\n', info.Modality);
fprintf('Series: %s\n', info.SeriesDescription);
fprintf('Spacing: %.2f x %.2f x %.2f mm\n', ...
    info.PixelSpacing(1), info.PixelSpacing(2), info.SliceThickness);

% Check volume size before loading
info = niftiinfo('large_scan.nii');
voxels = prod(info.ImageSize);
bytes_needed = voxels * info.BitsPerPixel / 8;
fprintf('Size: %.2f GB\n', bytes_needed / 1e9);

% Process slice-by-slice for huge volumes
for k = 1:V.NumTransverseSlices
    slice = extractSlice(V, k, 'transverse');
    % Process slice using IPT functions
    processed = imgaussfilt(slice, 1.5);  % IPT
    V = replaceSlice(V, processed, k, 'transverse');
end

% Load with MIT (preserves spatial info)
V = medicalVolume('scan.nii');

% Process with IPT (see matlab-image-processing-toolbox skill)
voxels = im2double(V.Voxels);
voxels = imgaussfilt3(voxels, 1.5);           % Gaussian smoothing
voxels = adapthisteq(voxels, 'NumTiles', [4 4 4]);  % CLAHE

% Segment with IPT
level = graythresh(voxels(:));
mask = voxels > level;
mask = imopen(mask, strel('sphere', 3));      % Morphological cleanup

% Create labeled volume with MIT (preserves spatial context)
labelVol = medicalVolume(uint8(mask), V.VolumeGeometry);
labelVol.Modality = 'SEG';
write(labelVol, 'segmentation.nii');

% From directory of DICOM files
V = medicalVolume('path/to/dicom_folder');

% From dicomCollection for multi-series
coll = dicomCollection('patient_folder', 'IncludeSubfolders', true);
disp(coll);  % Shows all series

% Load specific series
V = medicalVolume(coll, 'Rows', 2);  % 2nd series in collection

V = medicalVolume('brain.nii');

% Get middle slices in each plane
midT = round(V.NumTransverseSlices / 2);
midC = round(V.NumCoronalSlices / 2);
midS = round(V.NumSagittalSlices / 2);

transverse = extractSlice(V, midT, 'transverse');
coronal = extractSlice(V, midC, 'coronal');
sagittal = extractSlice(V, midS, 'sagittal');

% Display
figure;
subplot(1,3,1); imshow(transverse, []); title('Transverse');
subplot(1,3,2); imshow(coronal, []); title('Coronal');
subplot(1,3,3); imshow(sagittal, []); title('Sagittal');

fixed = medicalVolume('pre_contrast.nii');
moving = medicalVolume('post_contrast.nii');

% Rigid registration (fast, for same-modality)
[registered, tform] = imregmoment(moving, fixed);

% Affine registration with optimization
[optimizer, metric] = imregconfig('monomodal');
tform = imregtform(moving.Voxels, moving.VolumeGeometry, ...
                   fixed.Voxels, fixed.VolumeGeometry, ...
                   'affine', optimizer, metric);

% Apply transformation
registered = imwarp(moving.Voxels, tform, 'OutputView', ...
    imref3d(size(fixed.Voxels)));

V = medicalVolume('tumor_scan.nii');
mask = medicalVolume('tumor_mask.nii');

% Resample to isotropic (required for some features)
V_iso = resample(V, [1 1 1]);  % 1mm isotropic
mask_iso = resample(mask, [1 1 1]);

% Create radiomics object first (required pattern)
R = radiomics(V_iso.Voxels, mask_iso.Voxels > 0);

% Extract IBSI-compliant features as methods of radiomics object
intensity = intensityFeatures(R);
shape = shapeFeatures(R);
texture = textureFeatures(R);

% Combine into table
features = [intensity, shape, texture];

% Load model (requires support package)
model = medicalSegmentAnythingModel;

% Load image and extract embeddings
V = medicalVolume('liver_ct.nii');
slice = extractSlice(V, 50, 'transverse');
embeddings = extractEmbeddings(model, slice);

% Get image size (required parameter)
imageSize = size(slice);

% Segment with bounding box prompt (imageSize is required)
bbox = [100, 100, 200, 150];  % [x, y, width, height]
[mask, score] = segmentObjectsFromEmbeddings(model, embeddings, imageSize, ...
    BoundingBox=bbox);

% Combined workflow: MIT + IPT + Wavelet
V = medicalVolume('noisy_mri.nii');           % MIT: Load with spatial info

% Process each slice (IPT + Wavelet)
for k = 1:size(V.Voxels, 3)
    slice = im2double(V.Voxels(:,:,k));       % IPT: Convert type
    slice = wdenoise2(slice);                  % Wavelet: Denoise
    slice = adapthisteq(slice);                % IPT: Enhance
    V.Voxels(:,:,k) = slice;
end

write(V, 'enhanced_mri.nii');                  % MIT: Write with spatial info