diff --git a/2502_new/.gitignore b/2502_new/.gitignore new file mode 100644 index 0000000..36deac7 --- /dev/null +++ b/2502_new/.gitignore @@ -0,0 +1,22 @@ +# Python 缓存 +__pycache__/ +*.pyc + +# 虚拟环境 + #venv/ + +# IDE 设置（例如 PyCharm） +.idea/ + + # for ss2025 数据 +*.csv +*.xlsx +*.jpg +*.png +*.dcm +*.raw +*.mhd +*.nii +*.nii.gz +*.npz +*.attrs diff --git a/2502_new/2502-01.py b/2502_new/2502-01.py new file mode 100644 index 0000000..08e4e14 --- /dev/null +++ b/2502_new/2502-01.py @@ -0,0 +1,277 @@ +import argparse +import os +import numpy as np +import pydicom +from typing import List, Tuple +import SimpleITK as sitk + + +class DICOMProcessor: + def __init__(self): + self.dicom_files = [] + self.volume_data = None + self.metadata = {} + + def get_dicom_files(self, folder_path: str) -> List[str]: + """Get list of DICOM files in specified folder""" + print(f"Scanning folder: {folder_path}") + dicom_files = [] + + for file in os.listdir(folder_path): + file_path = os.path.join(folder_path, file) + if file.lower().endswith('.dcm') and os.path.isfile(file_path): + dicom_files.append(file_path) + print(f"Found DICOM file: {file}") + + if not dicom_files: + raise ValueError(f"No DICOM files (*.dcm) found in folder {folder_path}") + + print(f"Total {len(dicom_files)} DICOM files found") + return dicom_files + + def get_instance_number(self, dicom_file: str) -> int: + """Read DICOM file and get Instance Number""" + try: + ds = pydicom.dcmread(dicom_file, stop_before_pixels=True) + instance_number = int(getattr(ds, 'InstanceNumber', 0)) + return instance_number + except Exception as e: + print(f"Error reading file {dicom_file}: {e}") + return 0 + + def calculate_slice_spacing(self, ds1: pydicom.Dataset, ds2: pydicom.Dataset) -> float: + """Calculate Z-direction resolution (slice spacing)""" + + # Method 1: Use SliceLocation difference + slice_loc1 = getattr(ds1, 'SliceLocation', None) + slice_loc2 = getattr(ds2, 'SliceLocation', None) + + if slice_loc1 is not None and slice_loc2 is not None: + try: + spacing_from_slice_loc = abs(float(slice_loc2) - float(slice_loc1)) + print(f"Spacing calculated from SliceLocation: {spacing_from_slice_loc:.4f} mm") + return spacing_from_slice_loc + except (ValueError, TypeError): + pass + + # Method 2: Use Euclidean distance from ImagePositionPatient + pos1 = getattr(ds1, 'ImagePositionPatient', None) + pos2 = getattr(ds2, 'ImagePositionPatient', None) + + if pos1 is not None and pos2 is not None: + try: + pos1_array = np.array([float(pos1[0]), float(pos1[1]), float(pos1[2])]) + pos2_array = np.array([float(pos2[0]), float(pos2[1]), float(pos2[2])]) + spacing_from_position = np.linalg.norm(pos2_array - pos1_array) + print(f"Spacing calculated from ImagePositionPatient: {spacing_from_position:.4f} mm") + return spacing_from_position + except (ValueError, TypeError, IndexError): + pass + + print("Warning: Unable to calculate slice spacing, using default value 1.0 mm") + return 1.0 + + def extract_metadata(self, dicom_files: List[Tuple[str, int]]) -> dict: + """Extract metadata from DICOM files""" + if len(dicom_files) < 2: + print("Warning: Insufficient files to calculate slice spacing") + return {} + + file1, instance1 = dicom_files[0] + file2, instance2 = dicom_files[1] + + ds1 = pydicom.dcmread(file1) + ds2 = pydicom.dcmread(file2) + + metadata = {} + + # Get pixel spacing + metadata['pixel_spacing'] = getattr(ds1, 'PixelSpacing', [1.0, 1.0]) + print(f"Pixel spacing: {metadata['pixel_spacing']} mm") + + # Calculate slice spacing + metadata['slice_spacing'] = self.calculate_slice_spacing(ds1, ds2) + + # Get image orientation + metadata['image_orientation'] = getattr(ds1, 'ImageOrientationPatient', [1, 0, 0, 0, 1, 0]) + + # Get image position + metadata['image_position'] = getattr(ds1, 'ImagePositionPatient', [0, 0, 0]) + + # Get other metadata + metadata['rows'] = getattr(ds1, 'Rows', 0) + metadata['columns'] = getattr(ds1, 'Columns', 0) + metadata['bits_allocated'] = getattr(ds1, 'BitsAllocated', 16) + metadata['bits_stored'] = getattr(ds1, 'BitsStored', 12) + metadata['pixel_representation'] = getattr(ds1, 'PixelRepresentation', 0) + + return metadata + + def create_volume_data(self, folder_path: str) -> np.ndarray: + """Create 3D volume data""" + self.dicom_files = self.get_dicom_files(folder_path) + + print("Reading Instance Numbers...") + file_instance_pairs = [] + for file_path in self.dicom_files: + instance_number = self.get_instance_number(file_path) + file_instance_pairs.append((file_path, instance_number)) + + file_instance_pairs.sort(key=lambda x: x[1]) + print("Files sorted by Instance Number completed") + + self.metadata = self.extract_metadata(file_instance_pairs) + + print("Reading pixel data and creating volume data...") + slice_data = [] + for i, (file_path, instance_num) in enumerate(file_instance_pairs): + if i % 10 == 0 or i == len(file_instance_pairs) - 1: + print(f"Reading slice {i + 1}/{len(file_instance_pairs)}") + + ds = pydicom.dcmread(file_path) + pixel_array = ds.pixel_array + + if hasattr(ds, 'PixelRepresentation') and ds.PixelRepresentation == 1: + if hasattr(ds, 'RescaleSlope') and hasattr(ds, 'RescaleIntercept'): + pixel_array = pixel_array.astype(np.float32) * ds.RescaleSlope + ds.RescaleIntercept + + slice_data.append(pixel_array) + + volume_data = np.stack(slice_data, axis=0) + print(f"Volume data creation completed: {volume_data.shape}") + + return volume_data + + def save_as_mhd(self, volume_data: np.ndarray, output_filename: str): + """Save as MHD/RAW format""" + try: + image = sitk.GetImageFromArray(volume_data) + + if self.metadata: + pixel_spacing = self.metadata.get('pixel_spacing', [1.0, 1.0]) + slice_spacing = self.metadata.get('slice_spacing', 1.0) + image.SetSpacing([float(pixel_spacing[0]), float(pixel_spacing[1]), float(slice_spacing)]) + print(f"Setting voxel spacing: [{pixel_spacing[0]}, {pixel_spacing[1]}, {slice_spacing}]") + + direction = (1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0) + image.SetDirection(direction) + print("Setting direction matrix: [1, 0, 0, 0, 1, 0, 0, 0, -1]") + + sitk.WriteImage(image, output_filename) + print(f"MHD file saved: {output_filename}") + + self._enhance_mhd_file(output_filename) + + except Exception as e: + print(f"Error saving MHD file: {e}") + import traceback + traceback.print_exc() + + def _enhance_mhd_file(self, mhd_filename: str): + """Enhance MHD file by adding necessary metadata""" + try: + with open(mhd_filename, 'r') as f: + content = f.read() + + if 'TransformMatrix' not in content: + lines = content.split('\n') + new_lines = [] + + for line in lines: + new_lines.append(line) + if line.startswith('ElementSpacing'): + new_lines.append('TransformMatrix = 1 0 0 0 1 0 0 0 -1') + new_lines.append('Offset = 0 0 0') + new_lines.append('CenterOfRotation = 0 0 0') + new_lines.append('AnatomicalOrientation = RAS') + + with open(mhd_filename, 'w') as f: + f.write('\n'.join(new_lines)) + + print("MHD file metadata enhanced") + + except Exception as e: + print(f"Error enhancing MHD file: {e}") + + def save_as_nifti(self, volume_data: np.ndarray, output_filename: str): + """Save as NIfTI format""" + try: + image = sitk.GetImageFromArray(volume_data) + + if self.metadata: + pixel_spacing = self.metadata.get('pixel_spacing', [1.0, 1.0]) + slice_spacing = self.metadata.get('slice_spacing', 1.0) + image.SetSpacing([float(pixel_spacing[0]), float(pixel_spacing[1]), float(slice_spacing)]) + print(f"Setting voxel spacing: [{pixel_spacing[0]}, {pixel_spacing[1]}, {slice_spacing}]") + + sitk.WriteImage(image, output_filename) + print(f"NIfTI file saved: {output_filename}") + + except Exception as e: + print(f"Error saving NIfTI file: {e}") + import traceback + traceback.print_exc() + + +def main(): + parser = argparse.ArgumentParser(description='Convert DICOM files to 3D volume data') + parser.add_argument('input_folder', help='Input DICOM folder path') + parser.add_argument('output_file', help='Output file path (.mhd or .nii.gz)') + + args = parser.parse_args() + + # Check if input folder exists + if not os.path.isdir(args.input_folder): + print(f"Error: Folder {args.input_folder} does not exist") + return + + # Check output file extension + output_lower = args.output_file.lower() + if output_lower.endswith('.nii.gz'): + output_ext = '.nii.gz' + else: + output_ext = os.path.splitext(output_lower)[1] + + valid_extensions = ['.mhd', '.nii.gz'] + if output_ext not in valid_extensions: + print(f"Error: Unsupported output format. Supported formats: {valid_extensions}") + return + + try: + processor = DICOMProcessor() + + print("=" * 50) + print("Starting DICOM file processing...") + print("=" * 50) + + volume_data = processor.create_volume_data(args.input_folder) + + print("\n" + "=" * 50) + print("Volume data information:") + print("=" * 50) + print(f" Shape (Z, Y, X): {volume_data.shape}") + print(f" Data type: {volume_data.dtype}") + print(f" Value range: {np.min(volume_data):.2f} - {np.max(volume_data):.2f}") + + print("\n" + "=" * 50) + print("Saving file...") + print("=" * 50) + + # Save file based on extension + if output_ext == '.mhd': + processor.save_as_mhd(volume_data, args.output_file) + elif output_ext == '.nii.gz': + processor.save_as_nifti(volume_data, args.output_file) + + print("\n" + "=" * 50) + print("Conversion completed!") + print("=" * 50) + + except Exception as e: + print(f"\nError: Processing failed: {e}") + import traceback + traceback.print_exc() + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/2502_new/2502-02.py b/2502_new/2502-02.py new file mode 100644 index 0000000..89b4129 --- /dev/null +++ b/2502_new/2502-02.py @@ -0,0 +1,246 @@ +import SimpleITK as sitk +import numpy as np +from scipy import ndimage +from skimage import morphology +import argparse +import sys +import os + + +def read_nifti_file(file_path): + """Read NIfTI file and return image data and metadata""" + try: + image = sitk.ReadImage(file_path) + data = sitk.GetArrayFromImage(image) # Shape: (z, y, x) + return data, image + except Exception as e: + print(f"Error reading NIfTI file: {e}") + sys.exit(1) + + +def save_nifti_file(data, reference_image, output_path): + """Save NIfTI file while preserving original image geometry information""" + try: + output_image = sitk.GetImageFromArray(data) + output_image.CopyInformation(reference_image) + sitk.WriteImage(output_image, output_path) + print(f"Result saved to: {output_path}") + except Exception as e: + print(f"Error saving NIfTI file: {e}") + sys.exit(1) + + +def analyze_image_stats(data): + """Analyze image statistical information""" + print("=== Image Statistics ===") + print(f"Data type: {data.dtype}") + print(f"Image dimensions: {data.shape}") + print(f"Value range: [{np.min(data):.2f}, {np.max(data):.2f}]") + print(f"Mean: {np.mean(data):.2f}") + print(f"Standard deviation: {np.std(data):.2f}") + + # Display percentiles + percentiles = [0, 1, 5, 10, 25, 50, 75, 90, 95, 99, 100] + percentile_values = np.percentile(data, percentiles) + for p, val in zip(percentiles, percentile_values): + print(f"{p}% percentile: {val:.2f}") + + +def brain_ct_threshold(data): + """Threshold processing specifically for brain CT""" + print("=== Brain CT Threshold Processing ===") + + # Analyze data range + data_min, data_max = np.min(data), np.max(data) + print(f"Data range: [{data_min:.2f}, {data_max:.2f}]") + + # Typical HU value ranges for brain CT: + # - Air: -1000 HU + # - Fat: -100 to -50 HU + # - Water: 0 HU + # - CSF: 15 HU + # - Gray matter: 20-40 HU + # - White matter: 20-40 HU + # - Bone: 400-1000+ HU + + # Method 1: Use typical brain tissue HU value range + if data_min < -100: # Likely HU values + print("CT HU values detected, using brain tissue threshold") + # Typical HU value range for brain tissue (gray + white matter) + brain_tissue_min = 20 + brain_tissue_max = 80 + print(f"Brain tissue HU range: [{brain_tissue_min}, {brain_tissue_max}]") + + thresholded = np.where((data >= brain_tissue_min) & (data <= brain_tissue_max), 1, 0) + + # Method 2: If data is already normalized, use adaptive threshold + else: + print("Using adaptive threshold") + # Try multiple threshold strategies + threshold_strategies = [ + np.percentile(data, 60), # 60th percentile + np.percentile(data, 70), # 70th percentile + np.mean(data) + 0.5 * np.std(data), # Mean + 0.5 standard deviation + ] + + best_threshold = None + best_foreground_ratio = 0 + + for threshold in threshold_strategies: + temp_binary = (data >= threshold).astype(np.uint8) + foreground_ratio = np.sum(temp_binary) / temp_binary.size + + # Ideal brain region should occupy 20-60% of the image + if 0.2 <= foreground_ratio <= 0.6: + best_threshold = threshold + best_foreground_ratio = foreground_ratio + break + + if best_threshold is None: + # If no ideal threshold found, use 70th percentile + best_threshold = np.percentile(data, 70) + print(f"Using default threshold: {best_threshold:.2f} (70th percentile)") + else: + print(f"Selected best threshold: {best_threshold:.2f}, foreground ratio: {best_foreground_ratio:.2%}") + + thresholded = (data >= best_threshold).astype(np.uint8) + + print(f"Foreground voxels after thresholding: {np.sum(thresholded)}") + print(f"Foreground percentage after thresholding: {np.sum(thresholded) / thresholded.size * 100:.2f}%") + + return thresholded + + +def process_brain_slice(slice_data, slice_idx): + """Process single brain CT slice""" + # 1. Morphological processing - remove small noise + structure = np.ones((3, 3), dtype=np.uint8) + binary_slice = ndimage.binary_opening(slice_data, structure=structure).astype(np.uint8) + + # 2. Fill holes (especially in ventricular regions) + binary_slice = ndimage.binary_fill_holes(binary_slice).astype(np.uint8) + + # 3. Connected component analysis - using 8-connectivity + labeled_slice, num_features = ndimage.label(binary_slice, structure=np.ones((3, 3))) + + if num_features > 0: + # Calculate size of each component + component_sizes = [] + for i in range(1, num_features + 1): + component_sizes.append(np.sum(labeled_slice == i)) + + # Select the largest components (brain may have left and right hemispheres) + if len(component_sizes) >= 2: + # Take the top 2 largest components (left and right brain hemispheres) + largest_indices = np.argsort(component_sizes)[-2:] + binary_slice = np.zeros_like(slice_data) + for idx in largest_indices: + binary_slice = np.logical_or(binary_slice, labeled_slice == (idx + 1)) + elif len(component_sizes) >= 1: + # Only one large component + largest_component = np.argmax(component_sizes) + 1 + binary_slice = (labeled_slice == largest_component).astype(np.uint8) + + # 4. Final morphological cleanup + binary_slice = ndimage.binary_closing(binary_slice, structure=structure).astype(np.uint8) + binary_slice = ndimage.binary_fill_holes(binary_slice).astype(np.uint8) + + return binary_slice.astype(np.uint8) + + +def process_brain_volume_3d(binary_data, voxel_spacing): + """Brain 3D volume processing""" + # 1. 3D connected component analysis - using 6-connectivity + structure_3d = np.ones((3, 3, 3), dtype=np.uint8) + labeled_volume, num_features = ndimage.label(binary_data, structure=structure_3d) + + if num_features > 0: + # Calculate size of each 3D component + component_sizes = [] + for i in range(1, num_features + 1): + component_sizes.append(np.sum(labeled_volume == i)) + + # Select the largest component (brain) + if len(component_sizes) > 0: + largest_component = np.argmax(component_sizes) + 1 + binary_data = (labeled_volume == largest_component).astype(np.uint8) + + # 2. 3D morphological processing - considering anisotropy + z_ratio = voxel_spacing[2] / voxel_spacing[0] if voxel_spacing[0] > 0 else 1.0 + kernel_size_z = max(1, int(round(2 * z_ratio))) # Use smaller kernel + structure_3d_aniso = np.ones((kernel_size_z, 3, 3), dtype=np.uint8) + + # Gentle closing operation to fill small holes + binary_data = ndimage.binary_closing(binary_data, structure=structure_3d_aniso).astype(np.uint8) + + return binary_data + + +def extract_brain_region(input_file, output_file): + """Main function: Extract brain region""" + print(f"Input file: {input_file}") + + # Check if input file exists + if not os.path.exists(input_file): + print(f"Error: Input file does not exist: {input_file}") + sys.exit(1) + + # 1. Read image + data, image = read_nifti_file(input_file) + original_shape = data.shape + print(f"Image dimensions: {original_shape}") + + # Get voxel spacing information + voxel_spacing = image.GetSpacing() # (x, y, z) + print(f"Voxel spacing: {voxel_spacing}") + + # 2. Analyze image statistics + analyze_image_stats(data) + + # 3. Brain CT threshold processing + print("Performing brain CT threshold processing...") + binary_data = brain_ct_threshold(data) + + # Check threshold results + if np.sum(binary_data) == 0: + print("Warning: No foreground regions found after thresholding, trying more lenient threshold") + # Try more lenient threshold + if np.min(data) < -100: # HU values + binary_data = (data >= 0).astype(np.uint8) # All positive HU values + else: + binary_data = (data >= np.percentile(data, 40)).astype(np.uint8) # 40th percentile + print(f"Foreground voxels after lenient threshold: {np.sum(binary_data)}") + + # 4. Slice-by-slice processing + print("Performing slice-by-slice processing...") + processed_slices = np.zeros_like(binary_data, dtype=np.uint8) + + for z in range(binary_data.shape[0]): + if z % 20 == 0: # Show progress every 20 slices + print(f"Processing slice {z}/{binary_data.shape[0]}") + processed_slices[z] = process_brain_slice(binary_data[z], z) + + # 5. 3D volume processing + print("Performing 3D volume processing...") + final_mask = process_brain_volume_3d(processed_slices, voxel_spacing) + + # 6. Final statistics + print("=== Final Results ===") + print(f"Brain region voxel count: {np.sum(final_mask)}") + print(f"Brain region percentage: {np.sum(final_mask) / final_mask.size * 100:.2f}%") + + # 7. Save results + save_nifti_file(final_mask, image, output_file) + print("Brain extraction completed!") + + +# Direct execution using provided paths +if __name__ == "__main__": + input_path = "/Users/liulaolao/Desktop/R2/HeadCtSample_2022_volume 2.nii" + output_path = "/Users/liulaolao/Desktop/R2/brain_extraction_result.nii.gz" + + print("Starting brain CT image processing...") + print(f"Input file: {input_path}") + print(f"Output file: {output_path}") + + extract_brain_region(input_path, output_path) \ No newline at end of file diff --git a/2502_new/readme_01.md b/2502_new/readme_01.md new file mode 100644 index 0000000..30857ef --- /dev/null +++ b/2502_new/readme_01.md @@ -0,0 +1,95 @@ +# DICOM to 3D Volume Converter + +## Overview + +A Python tool for converting DICOM series into 3D volume data in MHD or NIfTI format. This tool automatically processes DICOM files, extracts metadata, calculates slice spacing, and generates properly formatted 3D medical images. + +## Features + +- ✅ Automatic DICOM file detection and sorting by Instance Number +- ✅ Slice spacing calculation using multiple methods +- ✅ Comprehensive metadata extraction +- ✅ Support for both MHD and NIfTI output formats +- ✅ Proper handling of DICOM rescaling (slope and intercept) +- ✅ Preservation of spatial information and orientation + +## Requirements + +### Python Dependencies +pydicom +numpy +SimpleITK + +text + +### Installation +```bash +pip install pydicom numpy SimpleITK +Usage + +Command Line Interface + +bash +# Convert to MHD format +python dicom_converter.py /path/to/dicom/folder output_volume.mhd + +# Convert to NIfTI format +python dicom_converter.py /path/to/dicom/folder output_volume.nii.gz +Parameters + +input_folder: Path to folder containing DICOM files (.dcm) +output_file: Output file path (.mhd or .nii.gz) +Processing Pipeline + +1. DICOM File Discovery + +Scans specified folder for .dcm files +Validates DICOM file integrity +Reports total number of files found +2. Instance Number Sorting + +Reads InstanceNumber from DICOM headers +Sorts files in correct slice order +Ensures proper volume reconstruction +3. Metadata Extraction + +Extracted Metadata: + +Pixel spacing (X, Y dimensions) +Slice spacing (Z dimension) +Image orientation and position +Rows and columns dimensions +Bits allocated and stored +Pixel representation +4. Slice Spacing Calculation + +Primary Methods: + +SliceLocation Difference: Uses SliceLocation DICOM tag +ImagePositionPatient: Calculates Euclidean distance between slices +Fallback: Defaults to 1.0 mm if calculation fails +5. Volume Data Creation + +Reads pixel arrays from all DICOM files +Applies rescaling (slope and intercept) if needed +Stacks slices into 3D volume array +Maintains original data type and values +6. File Export + +Supported Formats: + +MHD/RAW: MetaImage format with enhanced metadata +NIfTI: Standard neuroimaging format +Output Specifications + +MHD Format Features + +Proper voxel spacing in all three dimensions +Transform matrix for correct orientation +Anatomical orientation (RAS) +Offset and center of rotation information +NIfTI Format Features + +Correct spatial coordinates +Proper voxel dimensions +Compatible with most medical imaging software \ No newline at end of file diff --git a/2502_new/readme_02.md b/2502_new/readme_02.md new file mode 100644 index 0000000..3c706f7 --- /dev/null +++ b/2502_new/readme_02.md @@ -0,0 +1,192 @@ +# Brain CT/MR Extraction Tool + +## Overview + +A Python tool for automated brain region extraction from CT/MR images in NIfTI format. This tool processes medical imaging data to generate binary masks of brain regions, suitable for neuroimaging analysis and preprocessing pipelines. + +## Features + +- ✅ Support for both CT and MR imaging modalities +- ✅ Automatic thresholding and brain region extraction +- ✅ Morphological processing for noise removal and hole filling +- ✅ 3D connected component analysis +- ✅ Preservation of original image geometry +- ✅ Standard NIfTI format output + +## Requirements + +### Python Dependencies +``` +SimpleITK +numpy +scipy +scikit-image +``` + +### Installation +```bash +pip install SimpleITK numpy scipy scikit-image +``` + +## Usage + +### Method 1: Command Line Arguments +```bash +# Process CT images +python brain_extraction.py input_image.nii.gz output_mask.nii.gz --modality CT + +# Process MR images +python brain_extraction.py input_image.nii.gz output_mask.nii.gz --modality MR +``` + +### Method 2: Direct Code Execution +Modify the paths in the script and run directly: +```python +if __name__ == "__main__": + input_path = "/path/to/your/input.nii" + output_path = "/path/to/your/output.nii.gz" + modality = "CT" # or "MR" + + #extract_brain_region(input_path, output_path, modality) +``` + +## Processing Pipeline + +### 1. Image Reading & Analysis +- Read NIfTI format input images +- Analyze image statistics (value range, data type, percentiles) +- Extract voxel spacing information + +### 2. Threshold Processing +**CT Images**: +- Brain tissue HU range (20-80 HU) +- Adaptive thresholding strategies +- Multi-threshold validation + +**MR Images**: +- Otsu adaptive thresholding +- Slice-level threshold processing + +### 3. Morphological Processing +- 2D slice-wise opening operations for noise removal +- Hole filling (particularly for ventricular regions) +- 8-connected component analysis +- Preservation of major brain components + +### 4. 3D Volume Processing +- 6-connected 3D component analysis +- Anisotropic morphological processing +- Voxel spacing consideration +- Largest connected region extraction + +### 5. Output Generation +- Binary mask generation (brain=1, background=0) +- Preservation of original image geometry and spatial coordinates +- Output as uint8 NIfTI files + +## Output Specification + +Output files contain: +- **Brain region**: Pixel value = 1 +- **Background**: Pixel value = 0 +- **File format**: NIfTI (.nii or .nii.gz) +- **Data type**: np.uint8 + +## Examples + +### Process Brain CT Image +```bash +python brain_extraction.py \ + /Users/liulaolao/Desktop/R2/HeadCtSample_2022_volume_2.nii \ + /Users/liulaolao/Desktop/R2/brain_mask.nii.gz \ + --modality CT +``` + +### Result Verification +After processing, the program outputs: +- Image statistics +- Processing progress +- Final brain region percentage +- Output file path + +## Parameters + +### Command Line Arguments +- `input_file`: Input NIfTI file path +- `output_file`: Output NIfTI file path +- `--modality`: Imaging modality [CT|MR], default: CT + +### Threshold Parameters (CT Images) +- **Brain tissue HU range**: 20-80 HU +- **Window Width (WW)**: 50 +- **Window Level (WL)**: 400 +- **Adaptive thresholding**: Auto-adjusted based on image statistics + +## Troubleshooting + +### Common Issues + +1. **Completely Black Output** + - Verify input image is valid CT/MR data + - Confirm correct modality setting + - Check image statistics output + +2. **Insufficient Memory** + - Reduce number of simultaneously processed slices + - Use smaller structural elements + +3. **Long Processing Time** + - Consider downsampling for large images + - Adjust morphological processing parameters + +### Debug Information +The program provides detailed processing information: +- Image dimensions and data type +- Value ranges and statistics +- Thresholding results +- Processing progress at each stage +- Final extraction statistics + +## File Structure + +``` +brain_extraction.py # Main program file +README.md # Documentation +input_images/ # Input image directory +output_masks/ # Output mask directory +``` + +## Technical Details + +### Libraries Used +- **SimpleITK**: NIfTI file I/O operations +- **scipy.ndimage**: Morphological processing and connected component analysis +- **numpy**: Array operations and numerical computations +- **scikit-image**: Image processing algorithms + +### Algorithm Characteristics +- Automatic adaptation to different image contrasts +- Consideration of medical image anisotropy +- Robust noise handling +- Preservation of brain anatomical structure integrity + +## Version History + +- v1.0: Initial release with CT/MR brain extraction support +- v1.1: Optimized brain-specific thresholds, improved morphological processing + +## Important Notes + +1. Ensure input images are valid NIfTI format +2. Processing time may be longer for very large images +3. Recommend backing up original data before processing +4. Results should be verified using medical image viewers + +## Support + +For issues or suggestions, please check: +1. Input image format and integrity +2. Dependency package version compatibility +3. System memory availability +4. Output directory write permissions +``` \ No newline at end of file