OSCR

Accurate estimation of intravoxel incoherent motion parameters based on implicit neural representation.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches
  1. [1] § METHODS › Overview of the two‐stage IVIM‐INR framework ↔ train.py, lines 56–95 · score 0.76 · Xavier uniform initialization, hidden layers, ReLU, MLP, Weights, linear
  2. [2] § METHODS › Implementation details ↔ train.py, lines 450–490 · score 0.67 · cosine annealing learning, rate scheduling, Adam, GPU, Training, optimizer
  3. [3] § METHODS › Overview of the two‐stage IVIM‐INR framework ↔ train.py, lines 1–35 · score 0.60 · implicit neural representation, ReLU, IVIM INR, S2, global, S1
  4. [4] § METHODS › Overview of the two‐stage IVIM‐INR framework ↔ train.py, lines 56–95 · score 0.59 · layer MLP, hidden layers, ReLU, IVIM INR, SIREN, denoised

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,690 lines · 71 KB · no license · 4 matches

  1. #!/usr/bin/env python3
  2. """
  3. IVIM-INR: Two-stage implicit neural representation for IVIM parameter estimation.
  5. Stage 1 (S1): 4-layer ReLU MLP (512 hidden units per layer) that denoises each
  6. b-value image, taking as input spatial coordinates plus 3D patches
  7. from the other (non-target) b-value images.
  9. Stage 2 (S2): 4-layer SIREN (512 hidden units per layer) that maps spatial
  10. coordinates (plus a patch from the denoised b=0 image) to the
  11. IVIM parameters (D_p, D_t, F_p) and a global offset delta.
  13. Usage
  14. -----
  15. # Default: noisy SNR=50 brain phantom in ./data_1/SNR_50
  16. python ivim_fitting.py --data_dir data_1/SNR_50
  18. # Ground-truth (noise-free) fit
  19. python ivim_fitting.py --data_dir data_1 --no-use_noisy_data
  20. """
  22. import os
  23. import sys
  24. import time
  25. import json
  26. import argparse
  27. from typing import List, Tuple
  28. import numpy as np
  29. import nibabel as nib
  30. import torch
  31. import torch.nn as nn
  32. import torch.nn.functional as F
  33. from scipy import ndimage
  34. import pydicom
  35. from datetime import datetime
  37. # Configuration classes defined inline (previously from config_loader)
  38. # Default IVIM parameter ranges
  39. IVIM_PARAMETER_RANGES = {
  40. 'Dp': {'min': 0.001, 'max': 0.02}, # mm²/s
  41. 'Dt': {'min': 0.00001, 'max': 0.003}, # mm²/s
  42. 'Fp': {'min': 0.0, 'max': 0.5} # fraction
  43. }
  46. # ==================== Denoising Model Components ====================
  47. def fc_block(in_size, out_size, dropout, *args, **kwargs):
  48. """Fully connected block with ReLU and dropout."""
  49. return nn.Sequential(
  50. nn.Linear(in_size, out_size, *args, **kwargs),
  51. nn.ReLU(),
  52. nn.Dropout(dropout),
  53. )
  56. class MLPv1(nn.Module):
  57. """ReLU MLP used for DWI denoising (Stage 1 of IVIM-INR).
  59. num_layers counts TOTAL linear layers (input + hidden + output),
  60. matching the SirenNet convention. With the default num_layers=4 and
  61. hidden_size=512, this builds a 4-layer MLP where every hidden layer
  62. has 512 units (one input projection, two 512->512 hidden projections,
  63. and one output projection).
  64. """
  65. def __init__(self, input_size=3, hidden_size=512, output_size=1, dropout=0, num_layers=4):
  66. super(MLPv1, self).__init__()
  67. assert num_layers >= 2, "num_layers must be >= 2 (input + output)"
  68. self.input_size = input_size
  69. self.hidden_size = hidden_size
  70. self.output_size = output_size
  71. self.dropout = dropout
  72. self.num_layers = num_layers
  74. layers = [nn.Linear(input_size, hidden_size), nn.ReLU()]
  75. if dropout > 0:
  76. layers.append(nn.Dropout(dropout))
  77. for _ in range(num_layers - 2):
  78. layers.append(nn.Linear(hidden_size, hidden_size))
  79. layers.append(nn.ReLU())
  80. if dropout > 0:
  81. layers.append(nn.Dropout(dropout))
  82. layers.append(nn.Linear(hidden_size, output_size))
  84. # Xavier uniform initialization for all Linear layers
  85. for m in layers:
  86. if isinstance(m, nn.Linear):
  87. nn.init.xavier_uniform_(m.weight)
  88. if m.bias is not None:
  89. nn.init.zeros_(m.bias)
  91. self.net = nn.Sequential(*layers)
  93. def forward(self, x):
  94. x = x.view(-1, self.input_size)
  95. return self.net(x)
  98. class IntegratedPositionalEncoding(nn.Module):
  99. """Integrated Positional Encoding from MDINR."""
  100. def __init__(self, num_freqs: int, max_freq: float):
  101. super().__init__()
  102. self.freqs = torch.linspace(1.0, max_freq, steps=num_freqs)
  103. self.omega_max = max_freq
  105. def forward(self, x: torch.Tensor) -> torch.Tensor:
  106. """
  107. x: Tensor of shape [N, 3] (3D coordinates, normalized to [-1,1])
  108. Returns: Tensor of shape [N, 3 * F * 2] (concatenated sin/cos encoding)
  109. """
  110. N, D = x.shape
  111. device = x.device
  113. # Shape freqs to [1, 1, F]
  114. f = self.freqs.to(device).view(1, 1, -1)
  115. mu = x.view(N, D, 1)
  117. # Calculate sinc with d = π / ω_max
  118. d = torch.pi / self.omega_max
  119. sinc_arg = f * d
  121. sinc = torch.where(
  122. sinc_arg == 0,
  123. torch.ones_like(sinc_arg),
  124. torch.sin(sinc_arg) / sinc_arg
  125. )
  127. # Calculate sin and cos terms
  128. x_freqed = mu * f
  129. sin_term = torch.sin(x_freqed) * sinc
  130. cos_term = torch.cos(x_freqed) * sinc
  132. # Flatten and concatenate
  133. sin_term = sin_term.view(N, -1)
  134. cos_term = cos_term.view(N, -1)
  135. return torch.cat([sin_term, cos_term], dim=-1)
  138. # ==================== Scanner Type Detection ====================
  139. def detect_scanner_type(patient_id: str, date: str, dicom_base_dir: str) -> str:
  140. """
  141. Detect if the scanner is Unity (MR Linac) or MR Sim by checking DICOM header.
  143. Args:
  144. patient_id: Patient identifier
  145. date: Scan date
  146. dicom_base_dir: Base directory containing DICOM files
  148. Returns:
  149. 'unity' if Unity scanner, 'sim' if MR Sim scanner
  150. """
  151. # Find DWI directory
  152. patient_dir = os.path.join(dicom_base_dir, patient_id, date)
  154. if not os.path.exists(patient_dir):
  155. print(f" ⚠ Patient directory not found: {patient_dir}")
  156. return 'sim' # Default to sim
  158. # Look for DWI directories
  159. dwi_dirs = []
  160. for dirname in os.listdir(patient_dir):
  161. if 'DWI' in dirname.upper() and os.path.isdir(os.path.join(patient_dir, dirname)):
  162. dwi_dirs.append(os.path.join(patient_dir, dirname))
  164. if not dwi_dirs:
  165. print(f" ⚠ No DWI directories found in {patient_dir}")
  166. return 'sim' # Default to sim
  168. # Check first DICOM file in first DWI directory
  169. dwi_dir = dwi_dirs[0]
  170. dcm_files = [f for f in os.listdir(dwi_dir) if f.endswith('.dcm') or f.endswith('.DCM')]
  172. if not dcm_files:
  173. print(f" ⚠ No DICOM files found in {dwi_dir}")
  174. return 'sim' # Default to sim
  176. # Read first DICOM file
  177. dcm_path = os.path.join(dwi_dir, dcm_files[0])
  178. try:
  179. ds = pydicom.dcmread(dcm_path)
  181. # Check for (0008,1090) - ManufacturerModelName tag
  182. # Unity scanners have ManufacturerModelName = "Marlin"
  183. if hasattr(ds, 'ManufacturerModelName') or (0x0008, 0x1090) in ds:
  184. model_name = str(getattr(ds, 'ManufacturerModelName', '')).strip()
  185. print(f" ManufacturerModelName (0008,1090): {model_name}")
  187. if model_name.lower() == 'marlin':
  188. print(f" → Unity (MR Linac) detected")
  189. return 'unity'
  190. else:
  191. print(f" → MR Sim detected")
  192. return 'sim'
  193. else:
  194. print(f" No ManufacturerModelName tag (0008,1090) found - MR Sim detected")
  195. return 'sim'
  197. except Exception as e:
  198. print(f" ⚠ Error reading DICOM file {dcm_path}: {e}")
  199. return 'sim' # Default to sim
  202. # ==================== Noise Floor Detection ====================
  203. def detect_noise_floor_mask(b0_data: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
  204. """
  205. Detect noise floor regions in Unity data based on intensity thresholds.
  207. Args:
  208. b0_data: B0 image data
  210. Returns:
  211. background_mask: Pure background (intensity < 50)
  212. noise_floor_mask: Noise floor regions (50 <= intensity <= 300)
  213. tissue_mask: Signal regions (intensity > 300)
  214. """
  215. # Create masks based on thresholds
  216. background_mask = b0_data < 50
  217. noise_floor_mask = (b0_data >= 50) & (b0_data <= 300)
  218. tissue_mask = b0_data > 300
  220. # Get largest connected component of noise floor
  221. if np.any(noise_floor_mask):
  222. labeled, num_features = ndimage.label(noise_floor_mask)
  223. if num_features > 0:
  224. # Find largest component
  225. component_sizes = np.bincount(labeled.ravel())
  226. component_sizes[0] = 0 # Ignore background
  227. largest_component = component_sizes.argmax()
  229. # Keep only largest component as noise floor
  230. noise_floor_mask_clean = labeled == largest_component
  232. # Reassign small disconnected regions to tissue
  233. small_regions = noise_floor_mask & ~noise_floor_mask_clean
  234. tissue_mask = tissue_mask | small_regions
  235. noise_floor_mask = noise_floor_mask_clean
  237. return background_mask, noise_floor_mask, tissue_mask
  240. def rician_noise_correction(dwi_data: np.ndarray, noise_floor_mask: np.ndarray,
  241. method: str = "rayleigh_mean") -> Tuple[np.ndarray, float]:
  242. """
  243. Apply Rician noise floor correction to DWI data using noise floor regions.
  245. Args:
  246. dwi_data: DWI data array [X, Y, Z, B]
  247. noise_floor_mask: Binary mask indicating noise floor regions
  248. method: Correction method ('rayleigh_mean' or 'simple_subtraction')
  250. Returns:
  251. corrected_dwi: Corrected DWI data
  252. median_sigma: Estimated noise standard deviation
  253. """
  254. _, _, _, num_b = dwi_data.shape
  255. corrected_dwi = np.zeros_like(dwi_data)
  257. # Estimate noise from noise floor regions for each b-value
  258. sigma_values = []
  259. for b_idx in range(num_b):
  260. # Extract noise floor values
  261. noise_values = dwi_data[:, :, :, b_idx][noise_floor_mask]
  263. if len(noise_values) > 0:
  264. # For Rayleigh distribution: σ = mean / sqrt(π/2)
  265. mean_noise = np.mean(noise_values)
  266. sigma = mean_noise / np.sqrt(np.pi / 2)
  267. sigma_values.append(sigma)
  269. print(f" b={b_idx}: Noise floor mean={mean_noise:.3f}, σ={sigma:.3f}")
  271. if method == "rayleigh_mean":
  272. # Rician correction: S_corrected = sqrt(max(S^2 - 2σ^2, 0))
  273. S_squared = dwi_data[:, :, :, b_idx] ** 2
  274. S_corrected_squared = np.maximum(S_squared - 2 * sigma**2, 0)
  275. corrected_dwi[:, :, :, b_idx] = np.sqrt(S_corrected_squared)
  276. elif method == "simple_subtraction":
  277. # Simple subtraction method
  278. corrected_dwi[:, :, :, b_idx] = np.maximum(dwi_data[:, :, :, b_idx] - sigma, 0)
  279. else:
  280. print(f" b={b_idx}: No noise floor regions available, skipping correction")
  281. corrected_dwi[:, :, :, b_idx] = dwi_data[:, :, :, b_idx]
  282. sigma_values.append(0)
  284. # Use median sigma across b-values
  285. median_sigma = np.median(sigma_values) if sigma_values else 0
  286. print(f" Median σ across b-values: {median_sigma:.3f}")
  288. return corrected_dwi, median_sigma
  291. # ==================== Denoising Functions ====================
  292. def precompute_patch_features(noisy_dwi_padded: List[torch.Tensor], patch_coords: torch.Tensor,
  293. X_dim: int, Y_dim: int, Z_dim: int, num_b: int,
  294. device: torch.device, mask: np.ndarray = None) -> Tuple[torch.Tensor, np.ndarray]:
  295. """
  296. Precompute patch features for masked voxels and all b-values.
  298. Args:
  299. noisy_dwi_padded: List of padded DWI volumes for each b-value
  300. patch_coords: Patch coordinate offsets
  301. X_dim, Y_dim, Z_dim: Volume dimensions
  302. num_b: Number of b-values
  303. device: Device to use
  304. mask: Binary mask to select voxels (if None, use all voxels)
  306. Returns:
  307. all_patch_features: Tensor of shape [num_masked_voxels, num_b, patch_feature_size]
  308. voxel_coords: Array of voxel coordinates for masked voxels
  309. """
  310. print(" Precomputing patch features for masked voxels...")
  312. # Create coordinate grids
  313. all_indices = np.indices((X_dim, Y_dim, Z_dim))
  315. # Apply mask if provided
  316. if mask is not None:
  317. mask_indices = np.where(mask)
  318. voxel_coords = np.stack([mask_indices[0], mask_indices[1], mask_indices[2]], axis=1)
  319. print(f" Using mask: {len(voxel_coords):,} voxels out of {X_dim*Y_dim*Z_dim:,}")
  320. else:
  321. voxel_coords = np.stack([all_indices[0].ravel(),
  322. all_indices[1].ravel(),
  323. all_indices[2].ravel()], axis=1)
  325. num_voxels = len(voxel_coords)
  326. patch_feature_size = len(patch_coords)
  328. # Preallocate the feature tensor
  329. all_patch_features = torch.zeros((num_voxels, num_b, patch_feature_size),
  330. dtype=torch.float32, device=device)
  332. # Process in batches to manage memory
  333. batch_size = 5000
  335. for start_idx in range(0, num_voxels, batch_size):
  336. end_idx = min(start_idx + batch_size, num_voxels)
  337. batch_coords = voxel_coords[start_idx:end_idx]
  338. batch_coords_tensor = torch.tensor(batch_coords, dtype=torch.float32, device=device)
  340. # Extract patches for this batch
  341. # batch_coords are in original volume space, need to add padding offset
  342. # then add patch offsets (which already include padding)
  343. batch_coords_padded = batch_coords_tensor + 1 # Add padding offset to center coords
  344. coords_patch = batch_coords_padded[:, None, :] + patch_coords[None, :, :] - 1 # Subtract 1 because patch_coords already has padding
  345. x_idx = coords_patch[..., 0].reshape(-1).long()
  346. y_idx = coords_patch[..., 1].reshape(-1).long()
  347. z_idx = coords_patch[..., 2].reshape(-1).long()
  349. # Extract features for each b-value
  350. for b_idx in range(num_b):
  351. features = noisy_dwi_padded[b_idx][0, x_idx, y_idx, z_idx]
  352. features = features.reshape(len(batch_coords), -1)
  353. all_patch_features[start_idx:end_idx, b_idx, :] = features
  355. if (start_idx + batch_size) % 100000 == 0:
  356. print(f" Processed {end_idx}/{num_voxels} voxels")
  358. print(f" ✔ Precomputed patch features: shape {all_patch_features.shape}")
  359. print(f" Memory usage: {all_patch_features.element_size() * all_patch_features.nelement() / (1024**3):.2f} GB")
  361. return all_patch_features, voxel_coords
  364. def denoise_dwi_mdinr_optimized(dwi_data: np.ndarray, mask: np.ndarray, b_values: np.ndarray,
  365. args, device: torch.device, affine=None, save_intermediate=False,
  366. output_dir=None) -> np.ndarray:
  367. """
  368. Denoise DWI using MDINR INR model (Optimized version without DataLoader).
  370. Args:
  371. dwi_data: DWI data array [X, Y, Z, B]
  372. mask: Binary mask
  373. b_values: Array of b-values
  374. args: Command line arguments
  375. device: PyTorch device
  376. affine: Affine matrix for saving intermediate results (optional)
  377. save_intermediate: Whether to save each b-value after denoising (optional)
  378. output_dir: Directory for saving intermediate results (optional)
  380. Returns:
  381. denoised_dwi: Denoised DWI data
  382. """
  383. print(f"\n{'='*60}")
  384. print(f"DWI Denoising (MDINR INR - Optimized)")
  385. print(f"{'='*60}")
  387. X_dim, Y_dim, Z_dim, num_b = dwi_data.shape
  389. # Normalize DWI data by percentile
  390. percentiles = [np.percentile(dwi_data[:, :, :, b][mask], 99) for b in range(num_b)]
  391. noisy_dwi_norm = np.zeros_like(dwi_data)
  392. for b in range(num_b):
  393. noisy_dwi_norm[:, :, :, b] = dwi_data[:, :, :, b] / (percentiles[b] + 1e-6)
  395. # Setup positional encoding
  396. input_mapper = IntegratedPositionalEncoding(
  397. num_freqs=args.denoise_npe,
  398. max_freq=args.denoise_maxf
  399. ).to(device)
  401. # Convert to PyTorch tensors
  402. noisy_dwi_torch = [torch.from_numpy(noisy_dwi_norm[:, :, :, b]).float().unsqueeze(0).unsqueeze(0)
  403. for b in range(num_b)]
  405. # Pad volumes for patch extraction
  406. # For 3x3x3 patch, we need padding of 1 on each side
  407. pad_size = 1
  408. noisy_dwi_padded = [F.pad(img, (pad_size, pad_size, pad_size, pad_size, pad_size, pad_size),
  409. mode='constant', value=0) for img in noisy_dwi_torch]
  410. noisy_dwi_padded = [img.to(device)[0] for img in noisy_dwi_padded]
  412. # Create patch coordinates for 3x3x3 patch
  413. # Offsets should be [-1, 0, 1] in each dimension
  414. patch_offsets = [-1, 0, 1]
  415. patch_coords = np.array([(x, y, z) for x in patch_offsets for y in patch_offsets for z in patch_offsets])
  416. # Adjust for padding: add pad_size to make coordinates valid in padded volume
  417. patch_coords = patch_coords + pad_size
  418. patch_coords = torch.tensor(patch_coords, dtype=torch.float32).to(device)
  420. # Precompute all patch features (only for masked voxels)
  421. all_patch_features, masked_voxel_coords = precompute_patch_features(
  422. noisy_dwi_padded, patch_coords, X_dim, Y_dim, Z_dim, num_b, device, mask
  423. )
  425. # Process each b-value
  426. denoised_dwi = np.zeros_like(dwi_data)
  428. for b_idx in range(num_b):
  429. print(f"\nDenoising b-value {b_idx} (b={b_values[b_idx]})...")
  431. # Input size calculation
  432. # For 3x3x3 patch: 27 features per b-value, times (num_b - 1) excluding current b
  433. patch_feature_size = 27 * (num_b - 1)
  435. if args.denoise_use_coords:
  436. coord_encoding_size = args.denoise_npe * 2 * 3
  437. input_size = coord_encoding_size + patch_feature_size
  438. else:
  439. input_size = patch_feature_size
  441. # Create MLPv1 model
  442. model = MLPv1(
  443. input_size=input_size,
  444. output_size=1,
  445. hidden_size=args.denoise_hidden_size,
  446. num_layers=args.denoise_num_layers,
  447. dropout=args.denoise_dropout
  448. ).to(device)
  450. # Use masked voxel coordinates
  451. points = masked_voxel_coords.astype(np.float32)
  453. # Normalize coordinates
  454. points_normalized = np.empty_like(points, dtype=np.float32)
  455. points_normalized[:, 0] = 2.0 * points[:, 0] / (X_dim - 1) - 1.0
  456. points_normalized[:, 1] = 2.0 * points[:, 1] / (Y_dim - 1) - 1.0
  457. points_normalized[:, 2] = 2.0 * points[:, 2] / (Z_dim - 1) - 1.0
  459. # Move all data to GPU at once
  460. coords_normal = torch.tensor(points_normalized, dtype=torch.float32).to(device)
  461. coords_prior = torch.tensor(points, dtype=torch.float32).to(device)
  463. # Extract labels for masked voxels
  464. labels_masked = noisy_dwi_norm[:, :, :, b_idx][mask]
  465. labels = torch.tensor(labels_masked.reshape(-1, 1), dtype=torch.float32).to(device)
  467. if len(labels) == 0:
  468. continue
  470. # Get precomputed features for current b-value (excluding current b)
  471. current_b_features = torch.cat([
  472. all_patch_features[:, other_b, :]
  473. for other_b in range(num_b) if other_b != b_idx
  474. ], dim=1)
  476. # Training
  477. optimizer = torch.optim.Adam(model.parameters(), lr=args.denoise_lr)
  478. scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=args.denoise_epochs)
  479. criterion = nn.MSELoss()
  481. model.train()
  482. num_samples = len(coords_normal)
  483. batch_size = args.denoise_batch_size
  485. for epoch in range(args.denoise_epochs):
  486. epoch_loss = 0.0
  487. epoch_start_time = time.time()
  489. # Shuffle indices for this epoch
  490. indices = torch.randperm(num_samples, device=device)
  492. # Manual batching
  493. for start_idx in range(0, num_samples, batch_size):
  494. batch_start_time = time.time()
  496. end_idx = min(start_idx + batch_size, num_samples)
  497. batch_indices = indices[start_idx:end_idx]
  499. # Get batch data - already on GPU
  500. coords_normal_batch = coords_normal[batch_indices]
  501. # coords_prior_batch = coords_prior[batch_indices] # Not used in training
  502. label_batch = labels[batch_indices]
  503. patch_features_batch = current_b_features[batch_indices]
  505. # Build model input
  506. if args.denoise_use_coords:
  507. coords_encoded = input_mapper(coords_normal_batch)
  508. model_input = torch.cat([coords_encoded, patch_features_batch], dim=1)
  509. else:
  510. model_input = patch_features_batch
  512. # Forward pass
  513. pred = model(model_input)
  514. loss = criterion(pred, label_batch)
  516. # Backward pass
  517. optimizer.zero_grad()
  518. loss.backward()
  519. optimizer.step()
  521. epoch_loss += loss.item() * batch_indices.size(0)
  524. epoch_loss /= num_samples
  525. scheduler.step()
  527. if (epoch + 1) % 10 == 0:
  528. epoch_time = time.time() - epoch_start_time
  529. print(f" Epoch {epoch+1}/{args.denoise_epochs}, Loss: {epoch_loss:.6f}, Time: {epoch_time:.2f}s")
  531. # Inference
  532. model.eval()
  534. predictions = []
  535. with torch.no_grad():
  536. # Process in batches for inference
  537. for start_idx in range(0, num_samples, batch_size * 2): # Larger batch for inference
  538. end_idx = min(start_idx + batch_size * 2, num_samples)
  540. coords_norm_inf = coords_normal[start_idx:end_idx]
  541. patch_features_inf = current_b_features[start_idx:end_idx]
  543. # Build model input
  544. if args.denoise_use_coords:
  545. coords_encoded = input_mapper(coords_norm_inf)
  546. model_input = torch.cat([coords_encoded, patch_features_inf], dim=1)
  547. else:
  548. model_input = patch_features_inf
  550. pred = model(model_input).cpu().numpy()
  551. predictions.append(pred)
  553. # Reshape predictions - place them in the correct locations
  554. predictions = np.concatenate(predictions, axis=0)
  556. # Create full volume and place predictions at masked locations
  557. denoised_volume = np.zeros((X_dim, Y_dim, Z_dim), dtype=np.float32)
  558. for i, (x, y, z) in enumerate(masked_voxel_coords):
  559. denoised_volume[x, y, z] = predictions[i, 0] * percentiles[b_idx]
  561. denoised_dwi[:, :, :, b_idx] = denoised_volume
  563. # Save intermediate result if requested
  564. if save_intermediate and affine is not None and output_dir is not None:
  565. intermediate_path = os.path.join(output_dir, f'dwi_ivim_b{int(b_values[b_idx])}_denoised.nii.gz')
  566. if not os.path.exists(intermediate_path):
  567. intermediate_img = nib.Nifti1Image(denoised_volume, affine)
  568. nib.save(intermediate_img, intermediate_path)
  569. print(f" ✔ Saved intermediate denoised b={int(b_values[b_idx])} to: {intermediate_path}")
  570. else:
  571. print(f" ⏭ Skipping b={int(b_values[b_idx])} (already exists): {intermediate_path}")
  573. print("\n✔ Denoising completed")
  574. return denoised_dwi
  577. # ==================== Neural Network Components ====================
  578. class SirenLayer(nn.Module):
  579. """SIREN layer with sine activation."""
  580. def __init__(self, in_features: int, out_features: int,
  581. bias: bool = True, is_first: bool = False, omega_0: float = 30.0):
  582. super().__init__()
  583. self.omega_0 = omega_0
  584. self.is_first = is_first
  585. self.in_features = in_features
  586. self.linear = nn.Linear(in_features, out_features, bias=bias)
  587. self.init_weights()
  589. def init_weights(self):
  590. with torch.no_grad():
  591. if self.is_first:
  592. self.linear.weight.uniform_(-1 / self.in_features, 1 / self.in_features)
  593. else:
  594. bound = np.sqrt(6 / self.in_features) / self.omega_0
  595. self.linear.weight.uniform_(-bound, bound)
  597. def forward(self, x):
  598. return torch.sin(self.omega_0 * self.linear(x))
  601. class SirenNet(nn.Module):
  602. """Standard SIREN network (Sitzmann et al., NeurIPS 2020).
  604. num_layers counts TOTAL layers (first SIREN + hidden SIREN + final linear).
  605. """
  606. def __init__(self, input_size: int, output_size: int = 1,
  607. hidden_size: int = 512, num_layers: int = 4,
  608. first_omega_0: float = 30.0, hidden_omega_0: float = 30.0):
  609. super().__init__()
  610. assert num_layers >= 2, "num_layers must be >= 2"
  612. layers = []
  613. layers.append(SirenLayer(input_size, hidden_size,
  614. is_first=True, omega_0=first_omega_0))
  616. for _ in range(num_layers - 2):
  617. layers.append(SirenLayer(hidden_size, hidden_size,
  618. omega_0=hidden_omega_0))
  620. final_linear = nn.Linear(hidden_size, output_size)
  621. with torch.no_grad():
  622. bound = np.sqrt(6 / hidden_size) / hidden_omega_0
  623. final_linear.weight.uniform_(-bound, bound)
  625. layers.append(final_linear)
  626. self.net = nn.Sequential(*layers)
  628. def forward(self, x):
  629. return self.net(x)
  632. class IVIMNet(nn.Module):
  633. """Network for IVIM parameter estimation."""
  634. def __init__(self, backbone):
  635. super().__init__()
  636. self.backbone = backbone
  637. # Extract the output size from the last linear layer
  638. last_layer = backbone.net[-1]
  639. in_features = last_layer.in_features
  640. self.head = nn.Linear(in_features, 3)
  642. def forward(self, x):
  643. features = self.backbone.net[:-1](x) # All layers except last
  644. params = self.head(features)
  646. # Apply constraints to IVIM parameters using hardcoded ranges
  647. Dp = torch.sigmoid(params[:, 0:1]) * (IVIM_PARAMETER_RANGES['Dp']['max'] - IVIM_PARAMETER_RANGES['Dp']['min']) + IVIM_PARAMETER_RANGES['Dp']['min']
  648. Dt = torch.sigmoid(params[:, 1:2]) * (IVIM_PARAMETER_RANGES['Dt']['max'] - IVIM_PARAMETER_RANGES['Dt']['min']) + IVIM_PARAMETER_RANGES['Dt']['min']
  649. Fp = torch.sigmoid(params[:, 2:3]) * (IVIM_PARAMETER_RANGES['Fp']['max'] - IVIM_PARAMETER_RANGES['Fp']['min']) + IVIM_PARAMETER_RANGES['Fp']['min']
  651. return Dp, Dt, Fp
  654. class IVIMNetSoftConstraint(nn.Module):
  655. """Network for IVIM parameter estimation with soft b0 constraint."""
  656. def __init__(self, backbone, delta_min=-0.1, delta_max=0.1):
  657. super().__init__()
  658. self.backbone = backbone
  659. self.delta_min = delta_min
  660. self.delta_max = delta_max
  661. # Extract the output size from the last linear layer
  662. last_layer = backbone.net[-1]
  663. in_features = last_layer.in_features
  664. self.head = nn.Linear(in_features, 4) # Output 4 parameters: Dp, Dt, Fp, delta
  666. def forward(self, x):
  667. features = self.backbone.net[:-1](x) # All layers except last
  668. params = self.head(features)
  670. # Apply constraints to IVIM parameters using hardcoded ranges
  671. Dp = torch.sigmoid(params[:, 0:1]) * (IVIM_PARAMETER_RANGES['Dp']['max'] - IVIM_PARAMETER_RANGES['Dp']['min']) + IVIM_PARAMETER_RANGES['Dp']['min']
  672. Dt = torch.sigmoid(params[:, 1:2]) * (IVIM_PARAMETER_RANGES['Dt']['max'] - IVIM_PARAMETER_RANGES['Dt']['min']) + IVIM_PARAMETER_RANGES['Dt']['min']
  673. Fp = torch.sigmoid(params[:, 2:3]) * (IVIM_PARAMETER_RANGES['Fp']['max'] - IVIM_PARAMETER_RANGES['Fp']['min']) + IVIM_PARAMETER_RANGES['Fp']['min']
  674. delta = torch.tanh(params[:, 3:4]) * self.delta_max # tanh to get symmetric range [-delta_max, delta_max]
  676. return Dp, Dt, Fp, delta
  679. # ADCNet class removed - only IVIM fitting is performed
  682. # ==================== Signal Models ====================
  683. def compute_ivim_signal(Dp, Dt, Fp, b_vals):
  684. """Compute IVIM signal: Fp·exp(-b·Dp) + (1-Fp)·exp(-b·Dt)"""
  685. return Fp * torch.exp(-b_vals * Dp) + (1 - Fp) * torch.exp(-b_vals * Dt)
  688. def compute_ivim_signal_soft(Dp, Dt, Fp, delta, b_vals):
  689. """Compute IVIM signal with soft constraint: [Fp·exp(-b·Dp) + (1-Fp)·exp(-b·Dt)] + δ"""
  690. signal = Fp * torch.exp(-b_vals * Dp) + (1 - Fp) * torch.exp(-b_vals * Dt)
  691. return signal + delta
  694. # compute_adc_signal function removed - only IVIM fitting is performed
  697. # ==================== Parameter Fitting Functions ====================
  698. def fit_ivim_parameters(dwi_data: np.ndarray, mask: np.ndarray, b_values: np.ndarray,
  699. args, device: torch.device) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
  700. """
  701. Fit IVIM parameters using spatial SIREN network.
  703. Args:
  704. dwi_data: DWI data array [X, Y, Z, B]
  705. mask: Binary mask for fitting
  706. b_values: Array of b-values
  707. args: Command line arguments
  708. device: PyTorch device
  710. Returns:
  711. Dp_map, Dt_map, Fp_map: IVIM parameter maps
  712. """
  713. print(f"\n{'='*60}")
  714. print(f"IVIM Parameter Fitting")
  715. print(f"{'='*60}")
  717. # Filter b-values for IVIM based on hardcoded thresholds
  718. # Use b <= b_value_threshold and exclude b < min_b_value_threshold (except b=0)
  719. b_value_threshold = 1000 # Default threshold
  720. min_b_value_threshold = 10 # Default minimum threshold
  722. ivim_indices = np.where(
  723. (b_values <= b_value_threshold) &
  724. ((b_values >= min_b_value_threshold) | (b_values == 0))
  725. )[0]
  726. ivim_b_values = b_values[ivim_indices]
  727. ivim_dwi_data = dwi_data[:, :, :, ivim_indices]
  729. print(f" Using b-values <= {b_value_threshold} and >= {min_b_value_threshold} (except b=0) for IVIM fitting")
  730. print(f" Selected b-values: {ivim_b_values}")
  731. print(f" Number of b-values for IVIM: {len(ivim_b_values)}")
  733. X_dim, Y_dim, Z_dim, num_b_ivim = ivim_dwi_data.shape
  735. # Normalize by b0
  736. b0_expanded = ivim_dwi_data[:, :, :, 0:1]
  737. dwi_norm = np.clip(ivim_dwi_data / (b0_expanded + 1e-6), 0, 1)
  739. # Normalize B0 by its own maximum for use as additional input
  740. if args.ivim_use_b0_signal:
  741. b0_max = np.max(ivim_dwi_data[:, :, :, 0])
  742. b0_norm_by_max = ivim_dwi_data[:, :, :, 0] / (b0_max + 1e-6)
  743. print(f" B0 max value: {b0_max:.2f}, will use normalized B0 as additional input")
  745. # Prepare B0 for patch extraction if patch size > 1
  746. if args.ivim_b0_patch_size > 1:
  747. print(f" Using B0 patch size: {args.ivim_b0_patch_size}x{args.ivim_b0_patch_size}x{args.ivim_b0_patch_size}")
  748. # Convert to tensor and pad
  749. b0_tensor = torch.from_numpy(b0_norm_by_max).float().unsqueeze(0).unsqueeze(0)
  750. pad_size = args.ivim_b0_patch_size // 2
  751. b0_padded = F.pad(b0_tensor, (pad_size, pad_size, pad_size, pad_size, pad_size, pad_size),
  752. mode='constant', value=0)
  753. b0_padded = b0_padded.to(device)[0, 0] # Remove batch and channel dims
  755. # Create patch coordinates for B0
  756. patch_offsets = list(range(-pad_size, pad_size + 1))
  757. b0_patch_coords = np.array([(x, y, z) for x in patch_offsets for y in patch_offsets for z in patch_offsets])
  758. b0_patch_coords = b0_patch_coords + pad_size # Adjust for padding
  759. b0_patch_coords = torch.tensor(b0_patch_coords, dtype=torch.float32).to(device)
  760. b0_patch_size = len(b0_patch_coords)
  761. else:
  762. b0_patch_size = 1
  764. # Create coordinate grids
  765. X = np.arange(0, X_dim, 1)
  766. Y = np.arange(0, Y_dim, 1)
  767. Z = np.arange(0, Z_dim, 1)
  768. points = np.meshgrid(X, Y, Z, indexing='ij')
  769. points = np.stack(points).transpose(1, 2, 3, 0).reshape(-1, 3).astype(np.float32)
  771. # Normalize coordinates to [-1, 1]
  772. points_normalized = np.empty_like(points, dtype=np.float32)
  773. points_normalized[:, 0] = 2.0 * points[:, 0] / (X_dim - 1) - 1.0
  774. points_normalized[:, 1] = 2.0 * points[:, 1] / (Y_dim - 1) - 1.0
  775. points_normalized[:, 2] = 2.0 * points[:, 2] / (Z_dim - 1) - 1.0
  777. # Create IVIM model
  778. # Determine input size based on additional inputs
  779. input_size = 3 # Start with 3D coordinates
  781. if args.ivim_use_dwi_signals:
  782. input_size += num_b_ivim # Add normalized DWI signals
  783. print(f" Using DWI signals as additional input")
  785. if args.ivim_use_b0_signal:
  786. if args.ivim_b0_patch_size > 1:
  787. input_size += b0_patch_size # Add B0 patch features
  788. print(f" Using B0 patch features as additional input (size: {b0_patch_size})")
  789. else:
  790. input_size += 1 # Add single B0 value
  791. print(f" Using normalized B0 signal as additional input")
  793. print(f" Total input size: {input_size}")
  795. backbone = SirenNet(
  796. input_size=input_size,
  797. output_size=args.ivim_hidden_size,
  798. hidden_size=args.ivim_hidden_size,
  799. num_layers=args.ivim_num_layers,
  800. first_omega_0=args.ivim_first_omega_0,
  801. hidden_omega_0=args.ivim_hidden_omega_0,
  802. )
  804. model = IVIMNet(backbone).to(device)
  806. # Prepare data and move to GPU
  807. mask_flat = mask.reshape(-1)
  808. print(f" Mask voxels for fitting: {mask_flat.sum():,}")
  810. # Move masked data to GPU directly
  811. coords_tensor = torch.tensor(points_normalized[mask_flat], dtype=torch.float32).to(device)
  812. dwi_signals = torch.tensor(dwi_norm.reshape(-1, num_b_ivim)[mask_flat], dtype=torch.float32).to(device)
  814. if args.ivim_use_b0_signal:
  815. if args.ivim_b0_patch_size > 1:
  816. # Extract B0 patches for masked voxels
  817. print(f" Extracting B0 patches for {mask_flat.sum():,} masked voxels...")
  818. voxel_indices = np.where(mask_flat)[0]
  819. voxel_coords_full = points[voxel_indices]
  821. # Preallocate B0 patch features
  822. b0_patch_features = torch.zeros((len(voxel_coords_full), b0_patch_size), dtype=torch.float32, device=device)
  824. # Extract patches in batches
  825. batch_size = 5000
  826. for start_idx in range(0, len(voxel_coords_full), batch_size):
  827. end_idx = min(start_idx + batch_size, len(voxel_coords_full))
  828. batch_coords = voxel_coords_full[start_idx:end_idx]
  829. batch_coords_tensor = torch.tensor(batch_coords, dtype=torch.float32, device=device)
  831. # Extract B0 patches
  832. batch_coords_padded = batch_coords_tensor + pad_size
  833. coords_patch = batch_coords_padded[:, None, :] + b0_patch_coords[None, :, :] - pad_size
  834. x_idx = coords_patch[..., 0].reshape(-1).long()
  835. y_idx = coords_patch[..., 1].reshape(-1).long()
  836. z_idx = coords_patch[..., 2].reshape(-1).long()
  838. features = b0_padded[x_idx, y_idx, z_idx]
  839. features = features.reshape(len(batch_coords), -1)
  840. b0_patch_features[start_idx:end_idx, :] = features
  842. b0_signals = b0_patch_features
  843. else:
  844. b0_signals = torch.tensor(b0_norm_by_max.reshape(-1)[mask_flat], dtype=torch.float32).to(device).unsqueeze(1)
  846. # Exclude b0 from fitting (only use b>0)
  847. b_vals_t = torch.tensor(ivim_b_values[1:], device=device, dtype=torch.float32)
  849. # Training
  850. optimizer = torch.optim.Adam(model.parameters(), lr=args.ivim_lr)
  851. scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=args.ivim_epochs)
  852. criterion = nn.MSELoss()
  854. model.train()
  855. print(f"\nTraining for {args.ivim_epochs} epochs...")
  856. num_samples = len(coords_tensor)
  857. batch_size = args.ivim_batch_size
  859. for epoch in range(args.ivim_epochs):
  860. epoch_loss = 0.0
  861. epoch_start_time = time.time()
  863. # Shuffle indices for this epoch
  864. indices = torch.randperm(num_samples, device=device)
  866. # Manual batching
  867. for start_idx in range(0, num_samples, batch_size):
  868. end_idx = min(start_idx + batch_size, num_samples)
  869. batch_indices = indices[start_idx:end_idx]
  871. # Get batch data - already on GPU
  872. coords_batch = coords_tensor[batch_indices]
  873. signals_batch = dwi_signals[batch_indices]
  875. # Prepare model input
  876. model_input = coords_batch
  878. if args.ivim_use_dwi_signals:
  879. # Concatenate DWI signals
  880. model_input = torch.cat([model_input, signals_batch], dim=1)
  882. if args.ivim_use_b0_signal:
  883. # Concatenate B0 signal/patches
  884. b0_batch = b0_signals[batch_indices]
  885. model_input = torch.cat([model_input, b0_batch], dim=1)
  887. # Forward pass
  888. Dp, Dt, Fp = model(model_input)
  890. # Compute predicted signals (excluding b0)
  891. pred_signals = compute_ivim_signal(Dp, Dt, Fp, b_vals_t.unsqueeze(0))
  893. # Loss on normalized signals (excluding b0)
  894. loss = criterion(pred_signals, signals_batch[:, 1:])
  896. # Backward pass
  897. optimizer.zero_grad()
  898. loss.backward()
  899. torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
  900. optimizer.step()
  902. epoch_loss += loss.item() * batch_indices.size(0)
  904. epoch_loss /= num_samples
  905. scheduler.step()
  907. if (epoch + 1) % 10 == 0 or epoch == 0:
  908. epoch_time = time.time() - epoch_start_time
  909. print(f" Epoch {epoch+1}/{args.ivim_epochs}, Loss: {epoch_loss:.6f}, Time: {epoch_time:.2f}s")
  911. # Inference
  912. print("\nRunning inference...")
  913. model.eval()
  915. # Prepare full volume coordinates for inference
  916. coords_full = torch.tensor(points_normalized, dtype=torch.float32).to(device)
  917. nvox = coords_full.shape[0]
  918. step = 50000
  919. all_Dp, all_Dt, all_Fp = [], [], []
  921. with torch.no_grad():
  922. for i in range(0, nvox, step):
  923. coords_batch = coords_full[i:i+step]
  925. # Prepare model input for inference
  926. model_input = coords_batch
  928. if args.ivim_use_dwi_signals:
  929. # Get corresponding DWI signals for this batch
  930. dwi_signals_batch = torch.tensor(dwi_norm.reshape(-1, num_b_ivim)[i:i+step], dtype=torch.float32).to(device)
  931. model_input = torch.cat([model_input, dwi_signals_batch], dim=1)
  933. if args.ivim_use_b0_signal:
  934. if args.ivim_b0_patch_size > 1:
  935. # Extract B0 patches for inference batch
  936. batch_coords = points[i:i+step]
  937. batch_coords_tensor = torch.tensor(batch_coords, dtype=torch.float32, device=device)
  939. # Extract B0 patches
  940. batch_coords_padded = batch_coords_tensor + pad_size
  941. coords_patch = batch_coords_padded[:, None, :] + b0_patch_coords[None, :, :] - pad_size
  942. x_idx = coords_patch[..., 0].reshape(-1).long()
  943. y_idx = coords_patch[..., 1].reshape(-1).long()
  944. z_idx = coords_patch[..., 2].reshape(-1).long()
  946. features = b0_padded[x_idx, y_idx, z_idx]
  947. b0_signals_batch = features.reshape(coords_batch.shape[0], -1)
  948. else:
  949. # Get single B0 values
  950. b0_signals_batch = torch.tensor(b0_norm_by_max.reshape(-1)[i:i+step], dtype=torch.float32).to(device).unsqueeze(1)
  952. model_input = torch.cat([model_input, b0_signals_batch], dim=1)
  954. Dp, Dt, Fp = model(model_input)
  956. all_Dp.append(Dp.cpu().numpy())
  957. all_Dt.append(Dt.cpu().numpy())
  958. all_Fp.append(Fp.cpu().numpy())
  960. # Reshape to image dimensions and apply mask
  961. Dp_map = np.concatenate(all_Dp).reshape(X_dim, Y_dim, Z_dim) * mask
  962. Dt_map = np.concatenate(all_Dt).reshape(X_dim, Y_dim, Z_dim) * mask
  963. Fp_map = np.concatenate(all_Fp).reshape(X_dim, Y_dim, Z_dim) * mask
  965. print("✔ IVIM fitting completed")
  966. return Dp_map, Dt_map, Fp_map
  969. # fit_adc_parameters function removed - only IVIM fitting is performed
  972. def fit_ivim_parameters_soft(dwi_data: np.ndarray, mask: np.ndarray, b_values: np.ndarray,
  973. args, device: torch.device) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
  974. """
  975. Fit IVIM parameters using spatial SIREN network with soft b0 constraint.
  977. Args:
  978. dwi_data: DWI data array [X, Y, Z, B]
  979. mask: Binary mask for fitting
  980. b_values: Array of b-values
  981. args: Command line arguments
  982. device: PyTorch device
  983. fitting_config: Configuration object
  985. Returns:
  986. Dp_map, Dt_map, Fp_map, delta_map: IVIM parameter maps and delta offset map
  987. """
  988. print(f"\n{'='*60}")
  989. print(f"IVIM Parameter Fitting with Soft b0 Constraint")
  990. print(f"{'='*60}")
  992. # Filter b-values for IVIM based on hardcoded thresholds
  993. b_value_threshold = 1000 # Default threshold
  994. min_b_value_threshold = 10 # Default minimum threshold
  996. ivim_indices = np.where(
  997. (b_values <= b_value_threshold) &
  998. ((b_values >= min_b_value_threshold) | (b_values == 0))
  999. )[0]
  1000. ivim_b_values = b_values[ivim_indices]
  1001. ivim_dwi_data = dwi_data[:, :, :, ivim_indices]
  1003. print(f" Using b-values <= {b_value_threshold} and >= {min_b_value_threshold} (except b=0) for IVIM fitting")
  1004. print(f" Selected b-values: {ivim_b_values}")
  1005. print(f" Number of b-values for IVIM: {len(ivim_b_values)}")
  1007. X_dim, Y_dim, Z_dim, num_b_ivim = ivim_dwi_data.shape
  1009. # Normalize by b0 (soft constraint approach)
  1010. b0_expanded = ivim_dwi_data[:, :, :, 0:1]
  1011. dwi_norm = np.clip(ivim_dwi_data / (b0_expanded + 1e-6), 0, 2) # Allow some range above 1
  1013. print(f"\n DWI normalization for IVIM fitting:")
  1014. print(f" B0 range: [{b0_expanded[mask].min():.2f}, {b0_expanded[mask].max():.2f}]")
  1015. for i, b in enumerate(ivim_b_values):
  1016. print(f" b={b}: normalized range = [{dwi_norm[:,:,:,i][mask].min():.4f}, {dwi_norm[:,:,:,i][mask].max():.4f}]")
  1018. # Normalize B0 by its own maximum for use as additional input
  1019. if args.ivim_use_b0_signal:
  1020. b0_max = np.max(ivim_dwi_data[:, :, :, 0])
  1021. b0_norm_by_max = ivim_dwi_data[:, :, :, 0] / (b0_max + 1e-6)
  1022. print(f" B0 max value: {b0_max:.2f}, will use normalized B0 as additional input")
  1024. # Prepare B0 for patch extraction if patch size > 1
  1025. if args.ivim_b0_patch_size > 1:
  1026. print(f" Using B0 patch size: {args.ivim_b0_patch_size}x{args.ivim_b0_patch_size}x{args.ivim_b0_patch_size}")
  1027. # Convert to tensor and pad
  1028. b0_tensor = torch.from_numpy(b0_norm_by_max).float().unsqueeze(0).unsqueeze(0)
  1029. pad_size = args.ivim_b0_patch_size // 2
  1030. b0_padded = F.pad(b0_tensor, (pad_size, pad_size, pad_size, pad_size, pad_size, pad_size),
  1031. mode='constant', value=0)
  1032. b0_padded = b0_padded.to(device)[0, 0] # Remove batch and channel dims
  1034. # Create patch coordinates for B0
  1035. patch_offsets = list(range(-pad_size, pad_size + 1))
  1036. b0_patch_coords = np.array([(x, y, z) for x in patch_offsets for y in patch_offsets for z in patch_offsets])
  1037. b0_patch_coords = b0_patch_coords + pad_size # Adjust for padding
  1038. b0_patch_coords = torch.tensor(b0_patch_coords, dtype=torch.float32).to(device)
  1039. b0_patch_size = len(b0_patch_coords)
  1040. else:
  1041. b0_patch_size = 1
  1043. # Create coordinate grids
  1044. X = np.arange(0, X_dim, 1)
  1045. Y = np.arange(0, Y_dim, 1)
  1046. Z = np.arange(0, Z_dim, 1)
  1047. points = np.meshgrid(X, Y, Z, indexing='ij')
  1048. points = np.stack(points).transpose(1, 2, 3, 0).reshape(-1, 3).astype(np.float32)
  1050. # Normalize coordinates to [-1, 1]
  1051. points_normalized = np.empty_like(points, dtype=np.float32)
  1052. points_normalized[:, 0] = 2.0 * points[:, 0] / (X_dim - 1) - 1.0
  1053. points_normalized[:, 1] = 2.0 * points[:, 1] / (Y_dim - 1) - 1.0
  1054. points_normalized[:, 2] = 2.0 * points[:, 2] / (Z_dim - 1) - 1.0
  1056. # Create IVIM model with soft constraint
  1057. # Determine input size based on additional inputs
  1058. input_size = 3 # Start with 3D coordinates
  1060. if args.ivim_use_dwi_signals:
  1061. input_size += num_b_ivim # Add normalized DWI signals
  1062. print(f" Using DWI signals as additional input")
  1064. if args.ivim_use_b0_signal:
  1065. if args.ivim_b0_patch_size > 1:
  1066. input_size += b0_patch_size # Add B0 patch features
  1067. print(f" Using B0 patch features as additional input (size: {b0_patch_size})")
  1068. else:
  1069. input_size += 1 # Add single B0 value
  1070. print(f" Using normalized B0 signal as additional input")
  1072. print(f" Total input size: {input_size}")
  1074. backbone = SirenNet(
  1075. input_size=input_size,
  1076. output_size=args.ivim_hidden_size,
  1077. hidden_size=args.ivim_hidden_size,
  1078. num_layers=args.ivim_num_layers,
  1079. first_omega_0=args.ivim_first_omega_0,
  1080. hidden_omega_0=args.ivim_hidden_omega_0,
  1081. )
  1083. # Use soft constraint model
  1084. model = IVIMNetSoftConstraint(
  1085. backbone,
  1086. delta_min=args.delta_min,
  1087. delta_max=args.delta_max
  1088. ).to(device)
  1090. print(f"\n Delta (offset) constraint: [{args.delta_min:.2f}, {args.delta_max:.2f}]")
  1091. print(f" This allows the curve to deviate from exactly passing through (0,1)")
  1092. print(f" Delta regularization weight: {args.delta_regularization}")
  1094. # Prepare data and move to GPU
  1095. mask_flat = mask.reshape(-1)
  1096. print(f" Mask voxels for fitting: {mask_flat.sum():,}")
  1098. # Move masked data to GPU directly
  1099. coords_tensor = torch.tensor(points_normalized[mask_flat], dtype=torch.float32).to(device)
  1100. dwi_signals = torch.tensor(dwi_norm.reshape(-1, num_b_ivim)[mask_flat], dtype=torch.float32).to(device)
  1102. if args.ivim_use_b0_signal:
  1103. if args.ivim_b0_patch_size > 1:
  1104. # Extract B0 patches for masked voxels
  1105. print(f" Extracting B0 patches for {mask_flat.sum():,} masked voxels...")
  1106. voxel_indices = np.where(mask_flat)[0]
  1107. voxel_coords_full = points[voxel_indices]
  1109. # Preallocate B0 patch features
  1110. b0_patch_features = torch.zeros((len(voxel_coords_full), b0_patch_size), dtype=torch.float32, device=device)
  1112. # Extract patches in batches
  1113. batch_size = 5000
  1114. for start_idx in range(0, len(voxel_coords_full), batch_size):
  1115. end_idx = min(start_idx + batch_size, len(voxel_coords_full))
  1116. batch_coords = voxel_coords_full[start_idx:end_idx]
  1117. batch_coords_tensor = torch.tensor(batch_coords, dtype=torch.float32, device=device)
  1119. # Extract B0 patches
  1120. batch_coords_padded = batch_coords_tensor + pad_size
  1121. coords_patch = batch_coords_padded[:, None, :] + b0_patch_coords[None, :, :] - pad_size
  1122. x_idx = coords_patch[..., 0].reshape(-1).long()
  1123. y_idx = coords_patch[..., 1].reshape(-1).long()
  1124. z_idx = coords_patch[..., 2].reshape(-1).long()
  1126. features = b0_padded[x_idx, y_idx, z_idx]
  1127. features = features.reshape(len(batch_coords), -1)
  1128. b0_patch_features[start_idx:end_idx, :] = features
  1130. b0_signals = b0_patch_features
  1131. else:
  1132. b0_signals = torch.tensor(b0_norm_by_max.reshape(-1)[mask_flat], dtype=torch.float32).to(device).unsqueeze(1)
  1134. # Include ALL b-values for fitting with soft constraint (including b0)
  1135. b_vals_t = torch.tensor(ivim_b_values, device=device, dtype=torch.float32)
  1137. # Training
  1138. optimizer = torch.optim.Adam(model.parameters(), lr=args.ivim_lr)
  1139. scheduler = torch.optim.lr_scheduler.CosineAnnealingLR(optimizer, T_max=args.ivim_epochs)
  1141. model.train()
  1142. print(f"\nTraining for {args.ivim_epochs} epochs...")
  1143. num_samples = len(coords_tensor)
  1144. batch_size = args.ivim_batch_size
  1146. for epoch in range(args.ivim_epochs):
  1147. epoch_loss = 0.0
  1148. epoch_start_time = time.time()
  1150. # Shuffle indices for this epoch
  1151. indices = torch.randperm(num_samples, device=device)
  1153. # Manual batching
  1154. for start_idx in range(0, num_samples, batch_size):
  1155. end_idx = min(start_idx + batch_size, num_samples)
  1156. batch_indices = indices[start_idx:end_idx]
  1158. # Get batch data - already on GPU
  1159. coords_batch = coords_tensor[batch_indices]
  1160. signals_batch = dwi_signals[batch_indices]
  1162. # Prepare model input
  1163. model_input = coords_batch
  1165. if args.ivim_use_dwi_signals:
  1166. # Concatenate DWI signals
  1167. model_input = torch.cat([model_input, signals_batch], dim=1)
  1169. if args.ivim_use_b0_signal:
  1170. # Concatenate B0 signal/patches
  1171. b0_batch = b0_signals[batch_indices]
  1172. model_input = torch.cat([model_input, b0_batch], dim=1)
  1174. # Forward pass
  1175. Dp, Dt, Fp, delta = model(model_input)
  1177. # Compute predicted signals (including b0 with soft constraint)
  1178. pred_signals = compute_ivim_signal_soft(Dp, Dt, Fp, delta, b_vals_t.unsqueeze(0))
  1180. # Signal fitting loss
  1181. signal_loss = F.mse_loss(pred_signals, signals_batch)
  1183. # Delta regularization loss (L2 penalty on offset)
  1184. delta_reg_loss = args.delta_regularization * (delta**2).mean()
  1186. # Total loss
  1187. loss = signal_loss + delta_reg_loss
  1189. # Backward pass
  1190. optimizer.zero_grad()
  1191. loss.backward()
  1192. torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
  1193. optimizer.step()
  1195. epoch_loss += loss.item() * batch_indices.size(0)
  1197. epoch_loss /= num_samples
  1198. scheduler.step()
  1200. if (epoch + 1) % 10 == 0 or epoch == 0:
  1201. epoch_time = time.time() - epoch_start_time
  1202. print(f" Epoch {epoch+1}/{args.ivim_epochs}, Loss: {epoch_loss:.6f}, Time: {epoch_time:.2f}s")
  1204. # Inference
  1205. print("\nRunning inference...")
  1206. model.eval()
  1208. # Prepare full volume coordinates for inference
  1209. coords_full = torch.tensor(points_normalized, dtype=torch.float32).to(device)
  1210. nvox = coords_full.shape[0]
  1211. step = 50000
  1212. all_Dp, all_Dt, all_Fp, all_delta = [], [], [], []
  1214. with torch.no_grad():
  1215. for i in range(0, nvox, step):
  1216. coords_batch = coords_full[i:i+step]
  1218. # Prepare model input for inference
  1219. model_input = coords_batch
  1221. if args.ivim_use_dwi_signals:
  1222. # Get corresponding DWI signals for this batch
  1223. dwi_signals_batch = torch.tensor(dwi_norm.reshape(-1, num_b_ivim)[i:i+step], dtype=torch.float32).to(device)
  1224. model_input = torch.cat([model_input, dwi_signals_batch], dim=1)
  1226. if args.ivim_use_b0_signal:
  1227. if args.ivim_b0_patch_size > 1:
  1228. # Extract B0 patches for inference batch
  1229. batch_coords = points[i:i+step]
  1230. batch_coords_tensor = torch.tensor(batch_coords, dtype=torch.float32, device=device)
  1232. # Extract B0 patches
  1233. batch_coords_padded = batch_coords_tensor + pad_size
  1234. coords_patch = batch_coords_padded[:, None, :] + b0_patch_coords[None, :, :] - pad_size
  1235. x_idx = coords_patch[..., 0].reshape(-1).long()
  1236. y_idx = coords_patch[..., 1].reshape(-1).long()
  1237. z_idx = coords_patch[..., 2].reshape(-1).long()
  1239. features = b0_padded[x_idx, y_idx, z_idx]
  1240. b0_signals_batch = features.reshape(coords_batch.shape[0], -1)
  1241. else:
  1242. # Get single B0 values
  1243. b0_signals_batch = torch.tensor(b0_norm_by_max.reshape(-1)[i:i+step], dtype=torch.float32).to(device).unsqueeze(1)
  1245. model_input = torch.cat([model_input, b0_signals_batch], dim=1)
  1247. Dp, Dt, Fp, delta = model(model_input)
  1249. all_Dp.append(Dp.cpu().numpy())
  1250. all_Dt.append(Dt.cpu().numpy())
  1251. all_Fp.append(Fp.cpu().numpy())
  1252. all_delta.append(delta.cpu().numpy())
  1254. # Reshape to image dimensions and apply mask
  1255. Dp_map = np.concatenate(all_Dp).reshape(X_dim, Y_dim, Z_dim) * mask
  1256. Dt_map = np.concatenate(all_Dt).reshape(X_dim, Y_dim, Z_dim) * mask
  1257. Fp_map = np.concatenate(all_Fp).reshape(X_dim, Y_dim, Z_dim) * mask
  1258. delta_map = np.concatenate(all_delta).reshape(X_dim, Y_dim, Z_dim) * mask
  1260. print("\n IVIM parameter statistics:")
  1261. print(f" Dp range: [{Dp_map[mask].min():.6f}, {Dp_map[mask].max():.6f}] mm²/s")
  1262. print(f" Dt range: [{Dt_map[mask].min():.6f}, {Dt_map[mask].max():.6f}] mm²/s")
  1263. print(f" Fp range: [{Fp_map[mask].min():.4f}, {Fp_map[mask].max():.4f}]")
  1264. print(f" Delta range: [{delta_map[mask].min():.4f}, {delta_map[mask].max():.4f}]")
  1265. print(f" Mean delta: {delta_map[mask].mean():.4f} (mean absolute deviation: {np.abs(delta_map[mask]).mean():.4f})")
  1267. print("✔ IVIM fitting with soft b0 constraint completed")
  1268. return Dp_map, Dt_map, Fp_map, delta_map
  1271. # ==================== Main Processing Function ====================
  1272. def load_ivim_data(args):
  1273. """Load IVIM data from specified directory."""
  1274. print(f"\n{'#'*80}")
  1275. print(f"Loading IVIM data")
  1276. print(f"{'#'*80}")
  1278. # Build paths for IVIM data files
  1279. data_dir = args.data_dir if args.data_dir else args.input_dir
  1281. # Check if we should load noisy or ground truth data
  1282. if args.use_noisy_data:
  1283. data_subdir = os.path.join(data_dir, 'noisy')
  1284. print(f"Loading noisy IVIM data from: {data_subdir}")
  1285. else:
  1286. data_subdir = os.path.join(data_dir, 'ground_truth')
  1287. print(f"Loading ground truth IVIM data from: {data_subdir}")
  1289. # If subdirectory doesn't exist, fall back to main directory
  1290. if not os.path.exists(data_subdir):
  1291. data_subdir = data_dir
  1292. print(f"Subdirectory not found, using main directory: {data_dir}")
  1294. # Define b-values and corresponding files
  1295. b_values = [0, 50, 100, 200, 400, 600, 800, 1000]
  1296. dwi_files = {}
  1298. for b_val in b_values:
  1299. filename = f'dwi_ivim_b{b_val}.nii.gz'
  1300. filepath = os.path.join(data_subdir, filename)
  1301. if os.path.exists(filepath):
  1302. dwi_files[b_val] = filepath
  1303. else:
  1304. print(f" ⚠ File not found: {filepath}")
  1306. if not dwi_files:
  1307. print(f" ⚠ No IVIM files found in {data_subdir}")
  1308. return None, None, None
  1310. # Sort by b-value
  1311. b_values = sorted(dwi_files.keys())
  1312. print(f"Found {len(b_values)} b-values: {b_values}")
  1314. # Load DWI data
  1315. dwi_data = []
  1316. for b_val in b_values:
  1317. img = nib.load(dwi_files[b_val])
  1318. dwi_data.append(img.get_fdata().astype(np.float32))
  1320. affine = img.affine
  1321. dwi_data = np.stack(dwi_data, axis=-1) # Shape: [X, Y, Z, B]
  1322. b_values_array = np.array(b_values, dtype=np.float32)
  1324. print(f"DWI data shape: {dwi_data.shape}")
  1326. return dwi_data, b_values_array, affine
  1329. def process_ivim_fitting(args):
  1330. """Process IVIM parameter fitting."""
  1331. print(f"\n{'#'*80}")
  1332. print(f"IVIM parameter fitting")
  1333. print(f"{'#'*80}")
  1335. # This function has been replaced by process_ivim
  1337. # ADC fitting removed - only IVIM fitting is performed
  1340. def process_ivim(args):
  1341. """Process IVIM data (both loading and fitting)."""
  1342. # Validate output directory and files before processing
  1343. print(f"\nValidating output directory and file permissions...")
  1345. # Check if output directory exists, create if needed
  1346. try:
  1347. os.makedirs(args.output_dir, exist_ok=True)
  1348. except Exception as e:
  1349. sys.exit(f"Error: Cannot create output directory {args.output_dir}: {e}")
  1351. # Check write permissions
  1352. test_file = os.path.join(args.output_dir, '.test_write_permission')
  1353. try:
  1354. with open(test_file, 'w') as f:
  1355. f.write('test')
  1356. os.remove(test_file)
  1357. except Exception as e:
  1358. sys.exit(f"Error: Cannot write to output directory {args.output_dir}: {e}")
  1360. # Define expected output files
  1361. if args.use_soft_constraint:
  1362. output_files = ['Dp.nii.gz', 'Dt.nii.gz', 'Fp.nii.gz', 'delta.nii.gz', 'fitting_mask.nii.gz', 'parameter_config.json']
  1363. else:
  1364. output_files = ['Dp.nii.gz', 'Dt.nii.gz', 'Fp.nii.gz', 'fitting_mask.nii.gz', 'parameter_config.json']
  1366. # Check if output files already exist
  1367. existing_files = []
  1368. for filename in output_files:
  1369. filepath = os.path.join(args.output_dir, filename)
  1370. if os.path.exists(filepath):
  1371. existing_files.append(filename)
  1373. if existing_files and not args.force:
  1374. print(f"\n⚠ Warning: The following output files already exist:")
  1375. for f in existing_files:
  1376. print(f" - {f}")
  1377. return
  1379. # Check subdirectories if needed
  1380. if args.use_noisy_data and not args.skip_denoise:
  1381. denoised_dir = os.path.join(args.output_dir, 'denoised')
  1382. try:
  1383. os.makedirs(denoised_dir, exist_ok=True)
  1384. except Exception as e:
  1385. sys.exit(f"Error: Cannot create denoised directory {denoised_dir}: {e}")
  1387. if args.apply_noise_correction and args.use_noisy_data:
  1388. noise_corrected_dir = os.path.join(args.output_dir, 'noise_corrected')
  1389. try:
  1390. os.makedirs(noise_corrected_dir, exist_ok=True)
  1391. except Exception as e:
  1392. sys.exit(f"Error: Cannot create noise corrected directory {noise_corrected_dir}: {e}")
  1394. print("✔ Output directory validation completed\n")
  1396. # Load IVIM data
  1397. dwi_data, b_values_array, affine = load_ivim_data(args)
  1398. if dwi_data is None:
  1399. return
  1401. # Create fitting mask from b0
  1402. b0_data = dwi_data[:, :, :, 0]
  1404. # Load brain mask directly from data_1 directory (or the specific data directory being used)
  1405. # Use the input_dir or data_dir to construct mask path
  1406. base_data_dir = args.data_dir if args.data_dir else args.input_dir
  1407. # Extract the base directory (data_1, data_2, or data_3)
  1408. mask_path = os.path.join(os.path.dirname(base_data_dir) if 'SNR_' in base_data_dir or 'noisy' in base_data_dir or 'ground_truth' in base_data_dir else base_data_dir, 'brain_mask.nii.gz')
  1410. # Try multiple possible mask locations
  1411. possible_mask_paths = [
  1412. mask_path,
  1413. os.path.join(base_data_dir, 'brain_mask.nii.gz'),
  1414. os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data_1', 'brain_mask.nii.gz'),
  1415. ]
  1417. fitting_mask = None
  1418. for mp in possible_mask_paths:
  1419. if os.path.exists(mp):
  1420. mask_img = nib.load(mp)
  1421. fitting_mask = mask_img.get_fdata().astype(bool)
  1422. print(f"\nLoaded brain mask from: {mp}")
  1423. break
  1425. if fitting_mask is None:
  1426. # Fallback to creating mask from b0 data
  1427. fitting_mask = b0_data > args.background_threshold
  1428. print(f"\nWarning: Brain mask not found, using b0 threshold: {args.background_threshold}")
  1430. print(f"\nFitting mask coverage: {fitting_mask.sum() / fitting_mask.size * 100:.1f}%")
  1432. # Apply noise correction first if requested and using noisy data
  1433. if args.apply_noise_correction and args.use_noisy_data:
  1434. print("\nApplying noise floor correction to noisy data...")
  1436. # Detect noise floor regions from original b0
  1437. background_mask, noise_floor_mask, tissue_mask = detect_noise_floor_mask(b0_data)
  1439. # Apply Rician noise correction to original noisy data
  1440. dwi_data_corrected, estimated_sigma = rician_noise_correction(
  1441. dwi_data, noise_floor_mask, method=args.noise_correction_method
  1442. )
  1444. # Save noise corrected DWI data to separate folder
  1445. noise_corrected_dir = os.path.join(args.output_dir, 'noise_corrected')
  1446. os.makedirs(noise_corrected_dir, exist_ok=True)
  1448. # Save the noise-corrected data
  1449. for i, b_val in enumerate(b_values_array):
  1450. nib.save(nib.Nifti1Image(dwi_data_corrected[:, :, :, i].astype(np.float32), affine),
  1451. os.path.join(noise_corrected_dir, f'dwi_ivim_b{int(b_val)}_noise_corrected.nii.gz'))
  1453. # Use noise corrected data for subsequent processing
  1454. dwi_data = dwi_data_corrected
  1455. b0_data = dwi_data[:, :, :, 0]
  1456. print(f"\nNoise correction completed. Estimated sigma: {estimated_sigma:.4f}")
  1458. # Apply denoising after noise correction if requested and using noisy data
  1459. if args.use_noisy_data and not args.skip_denoise:
  1460. print("\nChecking for existing denoised data...")
  1461. device = torch.device(f'cuda:{args.gpu}' if torch.cuda.is_available() else 'cpu')
  1463. # Create denoised output directory
  1464. denoised_output_dir = os.path.join(args.output_dir, 'denoised')
  1465. os.makedirs(denoised_output_dir, exist_ok=True)
  1467. # Check if all denoised files exist
  1468. all_denoised_exist = True
  1469. denoised_data = []
  1470. for i, b_val in enumerate(b_values_array):
  1471. denoised_path = os.path.join(denoised_output_dir, f'dwi_ivim_b{int(b_val)}_denoised.nii.gz')
  1472. if os.path.exists(denoised_path):
  1473. print(f" ✔ Found existing denoised file for b={int(b_val)}")
  1474. # Load the existing denoised data
  1475. denoised_img = nib.load(denoised_path)
  1476. denoised_data.append(denoised_img.get_fdata().astype(np.float32))
  1477. else:
  1478. all_denoised_exist = False
  1479. break
  1481. if all_denoised_exist:
  1482. # All denoised files exist, use them
  1483. print(" ✔ All denoised files found, skipping denoising step")
  1484. dwi_data = np.stack(denoised_data, axis=-1)
  1485. else:
  1486. # Some denoised files missing, run denoising
  1487. print(" Some denoised files missing, running MDINR denoising...")
  1488. dwi_data = denoise_dwi_mdinr_optimized(dwi_data, fitting_mask, b_values_array, args, device,
  1489. affine=affine, save_intermediate=True,
  1490. output_dir=denoised_output_dir)
  1491. print("\nDenoising completed.")
  1493. # Update b0_data with denoised data
  1494. b0_data = dwi_data[:, :, :, 0]
  1496. elif args.use_noisy_data and args.skip_denoise:
  1497. print("\nUsing noisy data without denoising.")
  1498. else:
  1499. if not args.use_noisy_data:
  1500. print("\nUsing ground truth data (no noise correction or denoising needed).")
  1502. # Set device
  1503. device = torch.device(f'cuda:{args.gpu}' if torch.cuda.is_available() else 'cpu')
  1504. print(f"Using device: {device}")
  1506. # Configuration is now hardcoded in the script
  1507. # No need to load external config files
  1509. # Ensure output directory exists
  1510. os.makedirs(args.output_dir, exist_ok=True)
  1512. # Fit IVIM parameters
  1513. if args.use_soft_constraint:
  1514. Dp_map, Dt_map, Fp_map, delta_map = fit_ivim_parameters_soft(
  1515. dwi_data, fitting_mask, b_values_array, args, device
  1516. )
  1518. # Save IVIM parameters
  1519. nib.save(nib.Nifti1Image(Dp_map, affine), os.path.join(args.output_dir, 'Dp.nii.gz'))
  1520. nib.save(nib.Nifti1Image(Dt_map, affine), os.path.join(args.output_dir, 'Dt.nii.gz'))
  1521. nib.save(nib.Nifti1Image(Fp_map, affine), os.path.join(args.output_dir, 'Fp.nii.gz'))
  1522. nib.save(nib.Nifti1Image(delta_map, affine), os.path.join(args.output_dir, 'delta.nii.gz'))
  1523. print(f" ✔ Saved IVIM parameters with delta to: {args.output_dir}")
  1524. else:
  1525. Dp_map, Dt_map, Fp_map = fit_ivim_parameters(
  1526. dwi_data, fitting_mask, b_values_array, args, device
  1527. )
  1529. # Save IVIM parameters
  1530. nib.save(nib.Nifti1Image(Dp_map, affine), os.path.join(args.output_dir, 'Dp.nii.gz'))
  1531. nib.save(nib.Nifti1Image(Dt_map, affine), os.path.join(args.output_dir, 'Dt.nii.gz'))
  1532. nib.save(nib.Nifti1Image(Fp_map, affine), os.path.join(args.output_dir, 'Fp.nii.gz'))
  1533. print(f" ✔ Saved IVIM parameters to: {args.output_dir}")
  1535. # Save fitting mask
  1536. nib.save(nib.Nifti1Image(fitting_mask.astype(np.float32), affine),
  1537. os.path.join(args.output_dir, 'fitting_mask.nii.gz'))
  1539. # Save parameter configuration as JSON
  1540. config_info = {
  1541. 'IVIM_parameters': IVIM_PARAMETER_RANGES,
  1542. 'b_value_threshold': 1000,
  1543. 'min_b_value_threshold': 10,
  1544. 'soft_constraint': args.use_soft_constraint,
  1545. 'delta_min': args.delta_min if args.use_soft_constraint else None,
  1546. 'delta_max': args.delta_max if args.use_soft_constraint else None,
  1547. 'delta_regularization': args.delta_regularization if args.use_soft_constraint else None
  1548. }
  1549. config_path = os.path.join(args.output_dir, 'parameter_config.json')
  1550. with open(config_path, 'w') as f:
  1551. json.dump(config_info, f, indent=4)
  1552. print(f" ✔ Saved parameter configuration to: {config_path}")
  1554. print(f"\n✔ Completed IVIM fitting")
  1557. # process_patient function removed - direct IVIM processing only
  1559. # Removed - functionality merged into process_ivim
  1562. def main():
  1563. parser = argparse.ArgumentParser(description='IVIM Parameter Fitting Pipeline')
  1565. # Input/Output arguments
  1566. _this_dir = os.path.dirname(os.path.abspath(__file__))
  1567. parser.add_argument('--input_dir', type=str,
  1568. default=os.path.join(_this_dir, 'data_1'),
  1569. help='Input directory with IVIM data')
  1570. parser.add_argument('--data_dir', type=str,
  1571. default=None,
  1572. help='Data directory (overrides input_dir)')
  1573. parser.add_argument('--output_dir', type=str,
  1574. default=os.path.join(_this_dir, 'results', 'inr'),
  1575. help='Output directory for fitted parameters')
  1576. parser.add_argument('--denoised_dir', type=str,
  1577. default=None,
  1578. help='Output directory for denoised data (defaults to output_dir if not specified)')
  1579. parser.add_argument('--noise_corrected_dir', type=str,
  1580. default=None,
  1581. help='Output directory for noise corrected data (defaults to output_dir if not specified)')
  1583. # Processing options
  1584. parser.add_argument('--use_noisy_data', action='store_true', default=True,
  1585. help='Use noisy data instead of ground truth')
  1586. parser.add_argument('--apply_noise_correction', action='store_true', default=True,
  1587. help='Apply noise floor correction')
  1588. parser.add_argument('--noise_correction_method', choices=['rayleigh_mean', 'simple_subtraction'],
  1589. default='rayleigh_mean', help='Noise correction method')
  1590. parser.add_argument('--background_threshold', type=float, default=200,
  1591. help='Background threshold for creating masks')
  1592. parser.add_argument('--skip_denoise', action='store_true', help='Skip MDINR denoising step')
  1593. parser.add_argument('--force', action='store_true', help='Force reprocessing even if output files exist')
  1594. # Config parameter removed - configuration is now hardcoded
  1596. # Denoising parameters
  1597. parser.add_argument('--denoise_hidden_size', type=int, default=512)
  1598. parser.add_argument('--denoise_num_layers', type=int, default=4)
  1599. parser.add_argument('--denoise_dropout', type=float, default=0.0)
  1600. parser.add_argument('--denoise_npe', type=int, default=64)
  1601. parser.add_argument('--denoise_maxf', type=float, default=40)
  1602. parser.add_argument('--denoise_use_coords', type=int, default=0)
  1603. parser.add_argument('--denoise_epochs', type=int, default=800)
  1604. parser.add_argument('--denoise_batch_size', type=int, default=8000)
  1605. parser.add_argument('--denoise_lr', type=float, default=1e-4)
  1606. parser.add_argument('--denoise_patch_size', type=int, default=3, help='Patch size for denoising (will create NxNxN patches)')
  1608. # IVIM parameters
  1609. parser.add_argument('--ivim_hidden_size', type=int, default=512)
  1610. parser.add_argument('--ivim_num_layers', type=int, default=4)
  1611. parser.add_argument('--ivim_first_omega_0', type=float, default=30.0)
  1612. parser.add_argument('--ivim_hidden_omega_0', type=float, default=30.0)
  1613. parser.add_argument('--ivim_epochs', type=int, default=800)
  1614. parser.add_argument('--ivim_batch_size', type=int, default=6000)
  1615. parser.add_argument('--ivim_lr', type=float, default=5e-5)
  1616. parser.add_argument('--ivim_use_dwi_signals', action='store_true', default=False,
  1617. help='Use DWI signals as additional input for IVIM fitting')
  1618. parser.add_argument('--ivim_use_b0_signal', action='store_true', default=True,
  1619. help='Use normalized B0 signal as additional input for IVIM fitting')
  1620. parser.add_argument('--ivim_b0_patch_size', type=int, default=1,
  1621. help='Patch size for B0 input in IVIM fitting (default: 3, creates 3x3x3 patches)')
  1623. # Soft constraint parameters
  1624. parser.add_argument('--use_soft_constraint', action='store_true', default=True,
  1625. help='Use soft b0 constraint for IVIM fitting (allows curve to not pass exactly through (0,1))')
  1626. parser.add_argument('--delta_min', type=float, default=-0.01,
  1627. help='Minimum delta offset value for soft constraint')
  1628. parser.add_argument('--delta_max', type=float, default=0.01,
  1629. help='Maximum delta offset value for soft constraint')
  1630. parser.add_argument('--delta_regularization', type=float, default=0.01,
  1631. help='Regularization weight for delta offset (encourages delta close to 0)')
  1633. # Removed ADC parameters - only IVIM fitting is performed
  1635. # General parameters
  1636. parser.add_argument('--gpu', type=int, default=0, help='GPU device to use')
  1637. parser.add_argument('--seed', type=int, default=42, help='Random seed')
  1639. args = parser.parse_args()
  1641. # Set default directories if not provided
  1642. if args.denoised_dir is None:
  1643. args.denoised_dir = args.output_dir
  1644. if args.noise_corrected_dir is None:
  1645. args.noise_corrected_dir = args.output_dir
  1647. # Set random seeds
  1648. np.random.seed(args.seed)
  1649. torch.manual_seed(args.seed)
  1650. if torch.cuda.is_available():
  1651. torch.cuda.manual_seed(args.seed)
  1653. # Verify input directory
  1654. if not os.path.isdir(args.input_dir):
  1655. sys.exit(f'Error: Input directory not found: {args.input_dir}')
  1657. # Create output directory
  1658. os.makedirs(args.output_dir, exist_ok=True)
  1660. # No patient processing - direct IVIM data processing
  1662. print(f'='*80)
  1663. print('IVIM Parameter Fitting Pipeline')
  1664. print(f'='*80)
  1665. print(f'Input directory: {args.input_dir}')
  1666. print(f'Output directory: {args.output_dir}')
  1667. print(f'Use noisy data: {"Yes" if args.use_noisy_data else "No (Ground Truth)"}')
  1668. print(f'Denoising: {"Disabled" if args.skip_denoise else "Enabled"}')
  1669. print(f'Noise correction: {"Enabled" if args.apply_noise_correction else "Disabled"}')
  1670. print(f'Soft constraint: {"Enabled" if args.use_soft_constraint else "Disabled"}')
  1671. if args.use_soft_constraint:
  1672. print(f' Delta range: [{args.delta_min}, {args.delta_max}]')
  1673. print(f' Delta regularization: {args.delta_regularization}')
  1674. print(f'='*80)
  1676. # Process IVIM data
  1677. try:
  1678. process_ivim(args)
  1679. except Exception as e:
  1680. print(f'✖ Processing failed: {e}')
  1681. import traceback
  1682. traceback.print_exc()
  1684. print(f'\n{"="*80}')
  1685. print('Pipeline completed!')
  1686. print(f'{"="*80}')
  1689. if __name__ == '__main__':
  1690. main()

train.py at commit 9a92905, no license · at the source

Overview

Authors: Yunxiang Li1, Yen‐Peng Liao1, Yan Dai1, Jie Deng1, You Zhang1
  1. Department of Radiation Oncology, UT Southwestern Medical Center, Dallas, Texas, USA
Institutions: The University of Texas Southwestern Medical Center (United States)
Journal: Medical physics, volume 53, issue 8, article e70599
Dates: received 4 December 2025; accepted 1 July 2026; published online 28 July 2026; in print August 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1002/mp.70599 · PMID 42519880 · PMCID PMC13411821 · OpenAlex W7171561471
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: structural MRI / diffusion (modality), human (organism), other condition (population), systems (subfield)
Methods: Statistics, Preprocessing, Connectivity, Machine learning
Keywords: diffusion‐weighted imaging, implicit neural representation, intravoxel incoherent motion
MeSH: Diffusion Magnetic Resonance Imaging*, Image Processing, Computer-Assisted*, Brain, Humans, Imaging, Three-Dimensional, Motion, Movement, Phantoms, Imaging, Signal-To-Noise Ratio (* major topic)
Topic: MRI in cancer diagnosis (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: NCI NIH HHS (R01 CA258987, R01 CA240808, R01 CA280135); NIBIB NIH HHS (R01 EB034691); NIH HHS (R01 CA240808, R01 CA258987, R01 CA280135, R01 EB034691); National Institutes of Health (R01 CA240808, R01 CA280135, R01 EB034691, R01 CA258987)
Citations: not cited yet (Europe PMC); 52 references in the paper

Abstract

Background: Intravoxel incoherent motion (IVIM) diffusion‐weighted imaging has important value in treatment response monitoring. However, traditional voxel‐wise independent fitting methods are highly sensitive to noise and do not utilize spatial correlation, resulting in unstable parameter estimation. Although existing deep learning methods have shown improvements, they are still limited by local receptive fields.

Purpose: To address this, we propose a two‐stage IVIM parameter estimation framework based on Implicit Neural Representation (IVIM‐INR).

Methods: Our IVIM‐INR method achieves global spatial perception through coordinate encoding and enhances spatial context modeling by leveraging local 3D patch information from multi‐b‐value images. The first stage INR performs signal denoising, and the second stage INR accurately fits IVIM parameters.

Results: Evaluation on brain digital phantoms, AAPM breast IVIM‐dMRI Challenge data, and clinical Glioblastoma (GBM) patient data demonstrates significant advantages of the proposed method over existing techniques. In brain simulation data, when SNR = 50, the normalized mean absolute errors (NMAEs) in tumor regions were 0.16±0.10 for Dp, 0.02±0.01 for Dt, and 0.07±0.05 for Fp, all lower than comparison methods. In the tumor tissues of the 100 cases from the AAPM breast IVIM‐dMRI challenge dataset, Fp error was reduced by 58% compared to ConvNet. Intraclass correlation coefficient (ICC) analysis of real clinical data indicates that our method achieves the best performance with an ICC of Dt in normal tissues reaching 0.629.

Conclusions: By combining INR's continuous function modeling capability with spatial‐aware feature design, IVIM‐INR overcomes the inherent limitations of traditional methods under noisy conditions, providing a more reliable tool for clinical IVIM quantitative analysis.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

Kent0n-Li/IVIM-INR

License: none: the authors keep all their rights
State: the link answers, verified on 26 September 2026
Evidence: files inventoried
Commit: 9a92905c473df52d259d29f19aeec5e9f589e56d, 22 April 2026
Languages: Python (1)
Size: 2 files, 1 script
Software Heritage: not archived
Found in: the text, “Implementation details”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NiBabel (1 file), NumPy (1 file), pydicom (1 file), PyTorch (1 file), SciPy (1 file)
Availability: 1 check, the latest on 26 September 2026: the link answers
  • 26 September 2026: the link answers
2 files
At the source: github.com/Kent0n-Li/IVIM-INR

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 1 script, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 3, 28 September 2026

  • Publisher: n/a → Wiley

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 3 keywords, 9 MeSH terms, 4 funders, 45 references.

Cite

This paper

Li, Y., Liao, Y., Dai, Y., Deng, J., & Zhang, Y. (2026). Accurate estimation of intravoxel incoherent motion parameters based on implicit neural representation. Medical physics, 53(8), e70599. https://doi.org/10.1002/mp.70599

BibTeX

@article{li2026accurate,
author = {Li, Yunxiang and Liao, Yen‐Peng and Dai, Yan and Deng, Jie and Zhang, You},
title = {{Accurate estimation of intravoxel incoherent motion parameters based on implicit neural representation}},
journal = {Medical physics},
year = {2026},
month = aug,
volume = {53},
number = {8},
pages = {e70599},
publisher = {Wiley},
issn = {0094-2405},
doi = {10.1002/mp.70599},
url = {https://doi.org/10.1002/mp.70599},
pmid = {42519880},
pmcid = {PMC13411821}
}

RIS

TY - JOUR
AU - Li, Yunxiang
AU - Liao, Yen‐Peng
AU - Dai, Yan
AU - Deng, Jie
AU - Zhang, You
TI - Accurate estimation of intravoxel incoherent motion parameters based on implicit neural representation
T2 - Medical physics
J2 - Med Phys
PY - 2026
DA - 2026/08/01
VL - 53
IS - 8
SP - e70599
SN - 0094-2405
PB - Wiley
DO - 10.1002/mp.70599
UR - https://doi.org/10.1002/mp.70599
LA - en
ER -

CSL-JSON

{
"id": "10.1002/mp.70599",
"type": "article-journal",
"title": "Accurate estimation of intravoxel incoherent motion parameters based on implicit neural representation",
"container-title": "Medical physics",
"author": [
{
"family": "Li",
"given": "Yunxiang"
},
{
"family": "Liao",
"given": "Yen‐Peng"
},
{
"family": "Dai",
"given": "Yan"
},
{
"family": "Deng",
"given": "Jie"
},
{
"family": "Zhang",
"given": "You"
}
],
"container-title-short": "Med Phys",
"volume": "53",
"issue": "8",
"page": "e70599",
"DOI": "10.1002/mp.70599",
"PMID": "42519880",
"PMCID": "PMC13411821",
"ISSN": "0094-2405",
"publisher": "Wiley",
"URL": "https://doi.org/10.1002/mp.70599",
"language": "en",
"issued": {
"date-parts": [
[
2026,
8,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/nbm.70277 [code]
Hierarchical Bayesian Modelling Improves Microstructural Parameter Mapping in Diffusion and Exchange MRI Data.
Journal: NMR in biomedicine
In common: NiBabel, SciPy, NumPy, structural MRI / diffusion, 6 references
[2] doi:10.1038/s43856-026-01614-6 [code]
Simulation-based inference at the theoretical limit for fast, robust microstructural MRI with minimal diffusion data.
Journal: Communications medicine
In common: PyTorch, SciPy, NumPy, structural MRI / diffusion, 5 references
[3] doi:10.1002/mrm.70496 [code]
A Deep Nonlinear Subspace Modeling and Reconstruction for Diffusion-Weighted Imaging Using Denoising Auto-Encoder.
Journal: Magnetic resonance in medicine
In common: pydicom, NiBabel, PyTorch, 2 other tools, structural MRI / diffusion, 1 reference
[4] doi:10.1002/mrm.70405 [code]
DeepRelaxo: Fast Mono-Exponential Magnitude Brain R2* Mapping With Reduced Echoes Using Self-Supervised Deep Learning.
Journal: Magnetic resonance in medicine
In common: pydicom, NiBabel, PyTorch, 2 other tools, 1 reference
[5] doi:10.1080/07853890.2026.2685416 [code]
Pulmonary and cerebral damage in COVID-19 survivors: is there any association?
Journal: Annals of medicine
In common: pydicom, NiBabel, PyTorch, 2 other tools, structural MRI / diffusion, other condition
[6] doi:10.1002/alz.71530 [code]
Differential associations of plasma biomarkers with Alzheimer's disease and small vessel disease: A multimodal imaging study.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: NiBabel, PyTorch, SciPy, 1 other tool, structural MRI / diffusion, 2 references
[7] doi:10.3389/fnins.2026.1841093 [code]
Brain protein burden is related to intravoxel incoherent motion: PET-MR imaging study.
Journal: Frontiers in neuroscience
In common: NiBabel, SciPy, NumPy, structural MRI / diffusion, 2 references
[8] doi:10.1186/s13244-026-02296-3 [code]
A pre-trained foundation model framework for multiplanar MRI classification of extramural vascular invasion and mesorectal fascia invasion in rectal cancer.
Journal: Insights into imaging
In common: pydicom, NiBabel, PyTorch, 2 other tools, structural MRI / diffusion, other condition
[9] doi:10.1371/journal.pcbi.1014555 [code]
Body surface potential driven personalisation of electrophysiological digital twins in hypertrophic cardiomyopathy.
Journal: PLoS computational biology
In common: pydicom, NiBabel, PyTorch, 2 other tools, structural MRI / diffusion
[10] doi:10.3389/frai.2026.1771088 [code]
Few-shot deployment of pretrained MRI transformers in brain imaging tasks.
Journal: Frontiers in artificial intelligence
In common: pydicom, NiBabel, PyTorch, 2 other tools, structural MRI / diffusion

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.