3d_4d_viewer_pseudocode.txt

================================================================================
3D/4D RETINAL ANALYSIS SOFTWARE - DETAILED PSEUDOCODE
MVP Implementation Guide
================================================================================

OVERVIEW
--------
This pseudocode describes a minimal viable product for 3D/4D microscopy image
viewer similar to Imaris. The software loads volumetric image stacks (z-slices),
displays them as 3D renderings, supports time-series navigation, and provides
interactive camera controls.


================================================================================
1. DATA STRUCTURES
================================================================================

CLASS ImageVolume:
    PROPERTIES:
        width: INTEGER              // X dimension in pixels
        height: INTEGER             // Y dimension in pixels
        depth: INTEGER              // Z dimension (number of slices)
        voxelData: 3D_ARRAY[width][height][depth] of FLOAT  // Intensity values
        voxelSpacing: VECTOR3       // Physical spacing between voxels (x,y,z)
        bitDepth: INTEGER           // 8, 12, or 16 bit typically
        dataType: ENUM              // UINT8, UINT16, FLOAT32, etc.
    
    METHODS:
        GetVoxel(x, y, z) -> FLOAT
        SetVoxel(x, y, z, value)
        GetSlice(z) -> 2D_ARRAY


CLASS Channel:
    PROPERTIES:
        name: STRING                // e.g., "DAPI", "GFP", "RFP"
        volume: ImageVolume         // The 3D volume data
        color: RGB                  // Display color for this channel
        opacity: FLOAT              // 0.0 to 1.0
        intensityMin: FLOAT         // For contrast adjustment
        intensityMax: FLOAT         // For contrast adjustment
        visible: BOOLEAN            // Show/hide this channel
        transferFunction: ARRAY     // Mapping intensity to opacity
    
    METHODS:
        AdjustContrast(min, max)
        SetColor(r, g, b)
        SetOpacity(alpha)
        BuildTransferFunction()


CLASS TimePoint:
    PROPERTIES:
        timeIndex: INTEGER          // Frame number in sequence
        timestamp: FLOAT            // Actual time in seconds/minutes
        channels: LIST<Channel>     // All channels for this timepoint
    
    METHODS:
        GetChannel(index) -> Channel
        AddChannel(channel)


CLASS ImageDataset:
    PROPERTIES:
        timePoints: LIST<TimePoint> // All time frames
        numTimePoints: INTEGER
        numChannels: INTEGER
        metadata: DICTIONARY        // Microscope settings, date, etc.
        currentTimeIndex: INTEGER   // Currently displayed frame
    
    METHODS:
        GetCurrentTimePoint() -> TimePoint
        SetCurrentTimePoint(index)
        GetTimePoint(index) -> TimePoint


CLASS Camera:
    PROPERTIES:
        position: VECTOR3           // Camera position in world space
        target: VECTOR3             // Point camera is looking at
        up: VECTOR3                 // Up direction vector
        fieldOfView: FLOAT          // FOV in degrees
        nearPlane: FLOAT            // Near clipping plane
        farPlane: FLOAT             // Far clipping plane
        
        // Derived matrices
        viewMatrix: MATRIX4x4
        projectionMatrix: MATRIX4x4
    
    METHODS:
        UpdateViewMatrix()
        UpdateProjectionMatrix()
        Rotate(deltaX, deltaY)
        Pan(deltaX, deltaY)
        Zoom(deltaZ)
        LookAt(target, position, up)


CLASS VolumeRenderer:
    PROPERTIES:
        renderMode: ENUM            // RAY_CASTING, MAXIMUM_INTENSITY, etc.
        samplingRate: FLOAT         // Steps per voxel for ray marching
        quality: ENUM               // LOW, MEDIUM, HIGH
        enableShading: BOOLEAN
        lightDirection: VECTOR3
    
    METHODS:
        RenderVolume(dataset, camera) -> FRAMEBUFFER
        RayMarch(ray, volume) -> COLOR
        SampleVolume(position, volume) -> FLOAT
        ApplyTransferFunction(intensity, channel) -> COLOR


CLASS UserInterface:
    PROPERTIES:
        viewport: RECTANGLE         // 3D rendering area
        timeline: WIDGET            // Scrollbar for time navigation
        channelPanel: PANEL         // Channel visibility and color controls
        contrastPanel: PANEL        // Histogram and contrast sliders
        cameraControls: PANEL       // Rotation, zoom, reset buttons
    
    METHODS:
        Initialize()
        HandleMouseInput(event)
        HandleKeyboardInput(event)
        UpdateUI()


================================================================================
2. FILE IMPORT MODULE
================================================================================

FUNCTION LoadImageFile(filePath):
    """
    Load microscopy image files (TIFF, multi-page TIFF, proprietary formats)
    """
    
    // Detect file format
    fileExtension = GetFileExtension(filePath)
    
    IF fileExtension == ".tif" OR fileExtension == ".tiff":
        dataset = LoadTIFFStack(filePath)
    ELSE IF fileExtension == ".lsm":  // Zeiss format
        dataset = LoadLSM(filePath)
    ELSE IF fileExtension == ".czi":  // Zeiss CZI format
        dataset = LoadCZI(filePath)
    ELSE IF fileExtension == ".nd2":  // Nikon format
        dataset = LoadND2(filePath)
    ELSE:
        THROW Exception("Unsupported file format")
    
    RETURN dataset


FUNCTION LoadTIFFStack(filePath):
    """
    Load multi-page TIFF as 3D volume
    """
    
    // Open TIFF file
    tiffFile = OpenTIFF(filePath)
    
    // Read metadata from TIFF tags
    width = tiffFile.GetTag("ImageWidth")
    height = tiffFile.GetTag("ImageLength")
    numPages = tiffFile.GetPageCount()
    bitsPerSample = tiffFile.GetTag("BitsPerSample")
    samplesPerPixel = tiffFile.GetTag("SamplesPerPixel")
    
    // Determine if multi-channel or multi-z or both
    // Check for ImageJ/Fiji metadata
    imageDescription = tiffFile.GetTag("ImageDescription")
    metadata = ParseImageJMetadata(imageDescription)
    
    numChannels = metadata.channels OR 1
    numSlices = metadata.slices OR numPages
    numTimePoints = metadata.frames OR 1
    
    // Create dataset structure
    dataset = NEW ImageDataset()
    
    FOR t = 0 TO numTimePoints - 1:
        timePoint = NEW TimePoint()
        timePoint.timeIndex = t
        timePoint.timestamp = t * metadata.frameInterval
        
        FOR c = 0 TO numChannels - 1:
            channel = NEW Channel()
            channel.name = "Channel " + c
            volume = NEW ImageVolume()
            volume.width = width
            volume.height = height
            volume.depth = numSlices
            volume.voxelSpacing = VECTOR3(metadata.pixelWidth, 
                                          metadata.pixelHeight, 
                                          metadata.zSpacing)
            
            // Allocate memory for volume
            volume.voxelData = ALLOCATE_3D_ARRAY(width, height, numSlices)
            
            // Read each slice into the volume
            FOR z = 0 TO numSlices - 1:
                // Calculate page index: TIFF may be interleaved XYCZT
                pageIndex = CalculatePageIndex(c, z, t, numChannels, 
                                               numSlices, numTimePoints)
                
                tiffFile.SetCurrentPage(pageIndex)
                imageData = tiffFile.ReadImageData()
                
                // Copy slice into volume
                FOR y = 0 TO height - 1:
                    FOR x = 0 TO width - 1:
                        pixelValue = imageData[y][x]
                        volume.SetVoxel(x, y, z, pixelValue)
            
            channel.volume = volume
            channel.color = GetDefaultChannelColor(c)
            channel.opacity = 1.0
            channel.visible = TRUE
            channel.BuildTransferFunction()
            
            timePoint.AddChannel(channel)
        
        dataset.timePoints.ADD(timePoint)
    
    dataset.numTimePoints = numTimePoints
    dataset.numChannels = numChannels
    dataset.currentTimeIndex = 0
    
    tiffFile.Close()
    
    RETURN dataset


FUNCTION CalculatePageIndex(channel, slice, time, nChannels, nSlices, nTimes):
    """
    Calculate TIFF page index based on interleaving order
    Common orders: XYCZT, XYZCT, XYZTC
    """
    
    // Assume XYCZT order (most common)
    index = time * (nChannels * nSlices) + channel * nSlices + slice
    
    RETURN index


FUNCTION GetDefaultChannelColor(channelIndex):
    """
    Assign default colors to channels
    """
    
    colorPalette = [
        RGB(0, 0, 255),      // Blue (DAPI)
        RGB(0, 255, 0),      // Green (GFP)
        RGB(255, 0, 0),      // Red (RFP)
        RGB(255, 255, 0),    // Yellow
        RGB(255, 0, 255),    // Magenta
        RGB(0, 255, 255)     // Cyan
    ]
    
    RETURN colorPalette[channelIndex MOD colorPalette.LENGTH]


================================================================================
3. 3D VOLUME RENDERING
================================================================================

FUNCTION RenderVolume(dataset, camera, renderer):
    """
    Main rendering function using ray casting / ray marching
    """
    
    timePoint = dataset.GetCurrentTimePoint()
    
    // Get viewport dimensions
    viewportWidth = renderer.viewport.width
    viewportHeight = renderer.viewport.height
    
    // Create output framebuffer
    framebuffer = ALLOCATE_2D_ARRAY(viewportWidth, viewportHeight)
    
    // Update camera matrices
    camera.UpdateViewMatrix()
    camera.UpdateProjectionMatrix()
    
    // Calculate inverse view-projection matrix for ray generation
    invViewProj = INVERSE(camera.projectionMatrix * camera.viewMatrix)
    
    // For each pixel in viewport
    FOR y = 0 TO viewportHeight - 1:
        FOR x = 0 TO viewportWidth - 1:
            
            // Generate ray from camera through pixel
            ray = GenerateCameraRay(x, y, viewportWidth, viewportHeight, 
                                   camera, invViewProj)
            
            // Ray march through volume
            color = RayMarch(ray, timePoint, renderer)
            
            // Write to framebuffer
            framebuffer[y][x] = color
    
    RETURN framebuffer


FUNCTION GenerateCameraRay(pixelX, pixelY, screenWidth, screenHeight, 
                            camera, invViewProj):
    """
    Generate ray from camera through screen pixel
    """
    
    // Convert pixel coords to normalized device coordinates [-1, 1]
    ndcX = (2.0 * pixelX / screenWidth) - 1.0
    ndcY = 1.0 - (2.0 * pixelY / screenHeight)
    
    // Ray start at near plane
    rayStartNDC = VECTOR4(ndcX, ndcY, -1.0, 1.0)
    rayEndNDC = VECTOR4(ndcX, ndcY, 1.0, 1.0)
    
    // Transform to world space
    rayStartWorld = invViewProj * rayStartNDC
    rayEndWorld = invViewProj * rayEndNDC
    
    // Perspective divide
    rayStartWorld = rayStartWorld / rayStartWorld.w
    rayEndWorld = rayEndWorld / rayEndWorld.w
    
    // Create ray
    ray = NEW Ray()
    ray.origin = VECTOR3(rayStartWorld.x, rayStartWorld.y, rayStartWorld.z)
    ray.direction = NORMALIZE(VECTOR3(rayEndWorld.x - rayStartWorld.x,
                                      rayEndWorld.y - rayStartWorld.y,
                                      rayEndWorld.z - rayStartWorld.z))
    
    RETURN ray


FUNCTION RayMarch(ray, timePoint, renderer):
    """
    Ray marching through volume with compositing
    """
    
    // Calculate ray-box intersection with volume bounds
    volumeBounds = CalculateVolumeBounds(timePoint.channels[0].volume)
    
    tNear, tFar = IntersectRayBox(ray, volumeBounds)
    
    IF tNear < 0 OR tNear > tFar:
        // Ray misses volume
        RETURN RGB(0, 0, 0)
    
    // Initialize accumulation
    accumulatedColor = RGB(0, 0, 0)
    accumulatedOpacity = 0.0
    
    // Calculate step size based on sampling rate
    volume = timePoint.channels[0].volume
    maxDimension = MAX(volume.width, volume.height, volume.depth)
    stepSize = 1.0 / (maxDimension * renderer.samplingRate)
    
    // March along ray
    t = tNear
    WHILE t < tFar AND accumulatedOpacity < 0.95:
        
        // Calculate sample position
        samplePos = ray.origin + ray.direction * t
        
        // Sample all channels at this position
        FOR EACH channel IN timePoint.channels:
            IF channel.visible:
                // Trilinear interpolation
                intensity = SampleVolume(samplePos, channel.volume)
                
                // Apply transfer function
                sampleColor = channel.color * ApplyTransferFunction(intensity, 
                                                                     channel)
                sampleOpacity = ApplyTransferFunction(intensity, channel)
                
                // Modulate by channel opacity
                sampleOpacity = sampleOpacity * channel.opacity
                
                // Front-to-back compositing
                weight = sampleOpacity * (1.0 - accumulatedOpacity)
                accumulatedColor = accumulatedColor + sampleColor * weight
                accumulatedOpacity = accumulatedOpacity + weight
        
        // Advance along ray
        t = t + stepSize
    
    RETURN accumulatedColor


FUNCTION SampleVolume(position, volume):
    """
    Trilinear interpolation for smooth volume sampling
    """
    
    // Convert world position to voxel coordinates
    voxelX = position.x / volume.voxelSpacing.x
    voxelY = position.y / volume.voxelSpacing.y
    voxelZ = position.z / volume.voxelSpacing.z
    
    // Check bounds
    IF voxelX < 0 OR voxelX >= volume.width - 1 OR
       voxelY < 0 OR voxelY >= volume.height - 1 OR
       voxelZ < 0 OR voxelZ >= volume.depth - 1:
        RETURN 0.0
    
    // Get integer and fractional parts
    x0 = FLOOR(voxelX)
    y0 = FLOOR(voxelY)
    z0 = FLOOR(voxelZ)
    
    x1 = x0 + 1
    y1 = y0 + 1
    z1 = z0 + 1
    
    fx = voxelX - x0
    fy = voxelY - y0
    fz = voxelZ - z0
    
    // Sample 8 corners of voxel cube
    c000 = volume.GetVoxel(x0, y0, z0)
    c100 = volume.GetVoxel(x1, y0, z0)
    c010 = volume.GetVoxel(x0, y1, z0)
    c110 = volume.GetVoxel(x1, y1, z0)
    c001 = volume.GetVoxel(x0, y0, z1)
    c101 = volume.GetVoxel(x1, y0, z1)
    c011 = volume.GetVoxel(x0, y1, z1)
    c111 = volume.GetVoxel(x1, y1, z1)
    
    // Trilinear interpolation
    c00 = c000 * (1 - fx) + c100 * fx
    c01 = c001 * (1 - fx) + c101 * fx
    c10 = c010 * (1 - fx) + c110 * fx
    c11 = c011 * (1 - fx) + c111 * fx
    
    c0 = c00 * (1 - fy) + c10 * fy
    c1 = c01 * (1 - fy) + c11 * fy
    
    result = c0 * (1 - fz) + c1 * fz
    
    RETURN result


FUNCTION ApplyTransferFunction(intensity, channel):
    """
    Map intensity value to opacity using transfer function
    """
    
    // Normalize intensity based on contrast settings
    normalizedIntensity = (intensity - channel.intensityMin) / 
                          (channel.intensityMax - channel.intensityMin)
    
    // Clamp to [0, 1]
    normalizedIntensity = CLAMP(normalizedIntensity, 0.0, 1.0)
    
    // Simple linear transfer function for MVP
    // More complex: lookup table, piecewise functions, etc.
    opacity = normalizedIntensity
    
    RETURN opacity


FUNCTION IntersectRayBox(ray, bounds):
    """
    Calculate ray intersection with axis-aligned bounding box
    Returns near and far t values
    """
    
    // Slab method for AABB intersection
    tMin = (bounds.min.x - ray.origin.x) / ray.direction.x
    tMax = (bounds.max.x - ray.origin.x) / ray.direction.x
    
    IF tMin > tMax:
        SWAP(tMin, tMax)
    
    tyMin = (bounds.min.y - ray.origin.y) / ray.direction.y
    tyMax = (bounds.max.y - ray.origin.y) / ray.direction.y
    
    IF tyMin > tyMax:
        SWAP(tyMin, tyMax)
    
    IF tMin > tyMax OR tyMin > tMax:
        RETURN -1, -1  // No intersection
    
    tMin = MAX(tMin, tyMin)
    tMax = MIN(tMax, tyMax)
    
    tzMin = (bounds.min.z - ray.origin.z) / ray.direction.z
    tzMax = (bounds.max.z - ray.origin.z) / ray.direction.z
    
    IF tzMin > tzMax:
        SWAP(tzMin, tzMax)
    
    IF tMin > tzMax OR tzMin > tMax:
        RETURN -1, -1
    
    tMin = MAX(tMin, tzMin)
    tMax = MIN(tMax, tzMax)
    
    RETURN tMin, tMax


================================================================================
4. CAMERA CONTROLS
================================================================================

FUNCTION InitializeCamera(dataset):
    """
    Set up camera to view entire volume
    """
    
    camera = NEW Camera()
    
    // Get volume dimensions
    volume = dataset.GetTimePoint(0).GetChannel(0).volume
    
    // Calculate volume center
    centerX = (volume.width * volume.voxelSpacing.x) / 2.0
    centerY = (volume.height * volume.voxelSpacing.y) / 2.0
    centerZ = (volume.depth * volume.voxelSpacing.z) / 2.0
    
    volumeCenter = VECTOR3(centerX, centerY, centerZ)
    
    // Position camera at distance to view entire volume
    maxDimension = MAX(volume.width * volume.voxelSpacing.x,
                       volume.height * volume.voxelSpacing.y,
                       volume.depth * volume.voxelSpacing.z)
    
    cameraDistance = maxDimension * 2.0
    
    camera.position = volumeCenter + VECTOR3(0, 0, cameraDistance)
    camera.target = volumeCenter
    camera.up = VECTOR3(0, 1, 0)
    camera.fieldOfView = 45.0
    camera.nearPlane = 0.1
    camera.farPlane = cameraDistance * 10.0
    
    camera.UpdateViewMatrix()
    camera.UpdateProjectionMatrix()
    
    RETURN camera


FUNCTION HandleMouseRotate(camera, deltaX, deltaY, sensitivity):
    """
    Rotate camera around target using arcball rotation
    """
    
    // Convert mouse delta to angles
    angleX = deltaX * sensitivity
    angleY = deltaY * sensitivity
    
    // Calculate vector from target to camera
    toCamera = camera.position - camera.target
    distance = LENGTH(toCamera)
    
    // Create rotation around Y axis (horizontal rotation)
    rotationY = CreateRotationMatrix(camera.up, angleX)
    
    // Create right vector for vertical rotation
    forward = NORMALIZE(camera.target - camera.position)
    right = NORMALIZE(CROSS(forward, camera.up))
    
    // Create rotation around right axis (vertical rotation)
    rotationX = CreateRotationMatrix(right, angleY)
    
    // Apply rotations
    toCamera = rotationY * rotationX * toCamera
    
    // Update camera position
    camera.position = camera.target + toCamera
    
    // Update up vector (prevent camera from flipping)
    camera.up = NORMALIZE(rotationY * camera.up)
    
    camera.UpdateViewMatrix()


FUNCTION HandleMousePan(camera, deltaX, deltaY, sensitivity):
    """
    Pan camera parallel to view plane
    """
    
    // Calculate camera right and up vectors in world space
    forward = NORMALIZE(camera.target - camera.position)
    right = NORMALIZE(CROSS(forward, camera.up))
    up = NORMALIZE(CROSS(right, forward))
    
    // Calculate pan distance based on camera distance
    distance = LENGTH(camera.target - camera.position)
    panSpeed = distance * sensitivity
    
    // Calculate pan offset
    panOffset = right * (-deltaX * panSpeed) + up * (deltaY * panSpeed)
    
    // Move both camera and target
    camera.position = camera.position + panOffset
    camera.target = camera.target + panOffset
    
    camera.UpdateViewMatrix()


FUNCTION HandleMouseZoom(camera, delta, sensitivity):
    """
    Zoom camera toward/away from target
    """
    
    // Calculate vector from camera to target
    toTarget = camera.target - camera.position
    distance = LENGTH(toTarget)
    direction = toTarget / distance
    
    // Calculate zoom amount
    zoomAmount = delta * sensitivity * distance
    
    // Don't zoom too close or too far
    newDistance = distance + zoomAmount
    IF newDistance < 1.0:
        newDistance = 1.0
    IF newDistance > 10000.0:
        newDistance = 10000.0
    
    // Update camera position
    camera.position = camera.target - direction * newDistance
    
    camera.UpdateViewMatrix()


FUNCTION UpdateCameraMatrices(camera, viewportWidth, viewportHeight):
    """
    Update view and projection matrices
    """
    
    // Update view matrix (lookAt)
    camera.UpdateViewMatrix()
    
    // Update projection matrix (perspective)
    aspectRatio = viewportWidth / viewportHeight
    camera.projectionMatrix = CreatePerspectiveMatrix(
        camera.fieldOfView,
        aspectRatio,
        camera.nearPlane,
        camera.farPlane
    )


================================================================================
5. TIME NAVIGATION
================================================================================

FUNCTION SetupTimelineWidget(dataset, ui):
    """
    Initialize timeline scrubber for 4D navigation
    """
    
    timeline = ui.timeline
    timeline.minValue = 0
    timeline.maxValue = dataset.numTimePoints - 1
    timeline.currentValue = dataset.currentTimeIndex
    timeline.enabled = (dataset.numTimePoints > 1)
    
    // Display time information
    currentTime = dataset.GetCurrentTimePoint()
    timeline.label = "Time: " + currentTime.timestamp + " " + 
                     dataset.metadata.timeUnit


FUNCTION OnTimelineChanged(dataset, newTimeIndex, renderer):
    """
    Handle user changing time point
    """
    
    IF newTimeIndex < 0 OR newTimeIndex >= dataset.numTimePoints:
        RETURN
    
    // Update current time index
    dataset.currentTimeIndex = newTimeIndex
    
    // Trigger re-render with new time point
    renderer.needsUpdate = TRUE
    
    // Update UI
    UpdateTimelineDisplay(dataset)


FUNCTION PlayTimeSequence(dataset, renderer, frameRate):
    """
    Animate through time sequence
    """
    
    isPlaying = TRUE
    frameDelay = 1.0 / frameRate
    
    WHILE isPlaying:
        // Render current frame
        RenderFrame(dataset, renderer)
        
        // Advance to next time point
        dataset.currentTimeIndex = dataset.currentTimeIndex + 1
        
        // Loop back to beginning
        IF dataset.currentTimeIndex >= dataset.numTimePoints:
            dataset.currentTimeIndex = 0
        
        // Wait for frame delay
        SLEEP(frameDelay)
        
        // Check if user stopped playback
        IF userPressedStop:
            isPlaying = FALSE


================================================================================
6. MULTI-CHANNEL DISPLAY
================================================================================

FUNCTION SetupChannelPanel(dataset, ui):
    """
    Create UI controls for each channel
    """
    
    channelPanel = ui.channelPanel
    currentTime = dataset.GetCurrentTimePoint()
    
    FOR i = 0 TO currentTime.channels.LENGTH - 1:
        channel = currentTime.channels[i]
        
        // Create channel widget
        widget = NEW ChannelWidget()
        widget.label = channel.name
        widget.visibilityCheckbox = NEW Checkbox(channel.visible)
        widget.colorButton = NEW ColorPicker(channel.color)
        widget.opacitySlider = NEW Slider(0, 1, channel.opacity)
        
        // Add callbacks
        widget.visibilityCheckbox.OnChange = FUNCTION():
            channel.visible = widget.visibilityCheckbox.value
            renderer.needsUpdate = TRUE
        
        widget.colorButton.OnChange = FUNCTION():
            channel.color = widget.colorButton.selectedColor
            renderer.needsUpdate = TRUE
        
        widget.opacitySlider.OnChange = FUNCTION():
            channel.opacity = widget.opacitySlider.value
            renderer.needsUpdate = TRUE
        
        channelPanel.AddWidget(widget)


FUNCTION BlendChannels(channels, position):
    """
    Composite multiple channels at a sample position
    """
    
    finalColor = RGB(0, 0, 0)
    finalOpacity = 0.0
    
    FOR EACH channel IN channels:
        IF NOT channel.visible:
            CONTINUE
        
        // Sample this channel
        intensity = SampleVolume(position, channel.volume)
        
        // Apply transfer function
        opacity = ApplyTransferFunction(intensity, channel) * channel.opacity
        
        // Apply channel color
        channelColor = channel.color * intensity
        
        // Alpha blending
        weight = opacity * (1.0 - finalOpacity)
        finalColor = finalColor + channelColor * weight
        finalOpacity = finalOpacity + weight
        
        // Early exit if fully opaque
        IF finalOpacity >= 0.99:
            BREAK
    
    RETURN finalColor, finalOpacity


================================================================================
7. CONTRAST ADJUSTMENT
================================================================================

FUNCTION SetupContrastPanel(channel, ui):
    """
    Create histogram and contrast controls
    """
    
    contrastPanel = ui.contrastPanel
    
    // Calculate histogram of intensity values
    histogram = CalculateHistogram(channel.volume)
    
    // Display histogram
    histogramWidget = NEW HistogramWidget(histogram)
    contrastPanel.AddWidget(histogramWidget)
    
    // Min/max sliders
    minSlider = NEW Slider(0, 65535, channel.intensityMin)
    maxSlider = NEW Slider(0, 65535, channel.intensityMax)
    
    minSlider.OnChange = FUNCTION():
        channel.intensityMin = minSlider.value
        channel.BuildTransferFunction()
        renderer.needsUpdate = TRUE
    
    maxSlider.OnChange = FUNCTION():
        channel.intensityMax = maxSlider.value
        channel.BuildTransferFunction()
        renderer.needsUpdate = TRUE
    
    contrastPanel.AddWidget(minSlider)
    contrastPanel.AddWidget(maxSlider)
    
    // Auto contrast button
    autoButton = NEW Button("Auto Contrast")
    autoButton.OnClick = FUNCTION():
        AutoAdjustContrast(channel, histogram)
        minSlider.value = channel.intensityMin
        maxSlider.value = channel.intensityMax
        renderer.needsUpdate = TRUE
    
    contrastPanel.AddWidget(autoButton)


FUNCTION CalculateHistogram(volume):
    """
    Calculate intensity histogram for volume
    """
    
    numBins = 256
    histogram = ALLOCATE_ARRAY(numBins)
    FILL(histogram, 0)
    
    // Find min and max intensity
    minIntensity = INFINITY
    maxIntensity = -INFINITY
    
    FOR z = 0 TO volume.depth - 1:
        FOR y = 0 TO volume.height - 1:
            FOR x = 0 TO volume.width - 1:
                value = volume.GetVoxel(x, y, z)
                minIntensity = MIN(minIntensity, value)
                maxIntensity = MAX(maxIntensity, value)
    
    // Build histogram
    range = maxIntensity - minIntensity
    
    FOR z = 0 TO volume.depth - 1:
        FOR y = 0 TO volume.height - 1:
            FOR x = 0 TO volume.width - 1:
                value = volume.GetVoxel(x, y, z)
                
                // Map to bin
                normalized = (value - minIntensity) / range
                bin = FLOOR(normalized * (numBins - 1))
                histogram[bin] = histogram[bin] + 1
    
    RETURN histogram


FUNCTION AutoAdjustContrast(channel, histogram):
    """
    Auto-adjust contrast based on histogram (e.g., 1% and 99% percentiles)
    """
    
    volume = channel.volume
    totalVoxels = volume.width * volume.height * volume.depth
    
    // Find 1st and 99th percentile
    lowPercentile = 0.01
    highPercentile = 0.99
    
    lowThreshold = totalVoxels * lowPercentile
    highThreshold = totalVoxels * highPercentile
    
    cumulative = 0
    minBin = 0
    maxBin = histogram.LENGTH - 1
    
    FOR i = 0 TO histogram.LENGTH - 1:
        cumulative = cumulative + histogram[i]
        IF cumulative >= lowThreshold AND minBin == 0:
            minBin = i
        IF cumulative >= highThreshold:
            maxBin = i
            BREAK
    
    // Map bins back to intensity values
    // This requires knowing the original min/max used for histogram
    // For simplicity, assume 16-bit range
    channel.intensityMin = (minBin / histogram.LENGTH) * 65535
    channel.intensityMax = (maxBin / histogram.LENGTH) * 65535
    
    channel.BuildTransferFunction()


================================================================================
8. MAIN APPLICATION LOOP
================================================================================

FUNCTION Main():
    """
    Main application entry point
    """
    
    // Initialize graphics system
    InitializeGraphics()
    
    // Create main window
    window = CreateWindow("3D/4D Retinal Analysis Software", 1280, 720)
    
    // Initialize user interface
    ui = NEW UserInterface()
    ui.Initialize(window)
    
    // Create renderer
    renderer = NEW VolumeRenderer()
    renderer.renderMode = RAY_CASTING
    renderer.samplingRate = 2.0
    renderer.quality = MEDIUM
    
    // Dataset starts as null
    dataset = NULL
    camera = NULL
    
    // Variables for interaction
    mouseDown = FALSE
    lastMouseX = 0
    lastMouseY = 0
    
    // Main event loop
    WHILE window.isOpen:
        
        // Process events
        FOR EACH event IN GetEvents():
            
            IF event.type == FILE_OPEN:
                // User selected a file to open
                TRY:
                    dataset = LoadImageFile(event.filePath)
                    camera = InitializeCamera(dataset)
                    SetupTimelineWidget(dataset, ui)
                    SetupChannelPanel(dataset, ui)
                    renderer.needsUpdate = TRUE
                CATCH error:
                    ShowErrorDialog("Failed to load file: " + error.message)
            
            ELSE IF event.type == MOUSE_DOWN:
                mouseDown = TRUE
                lastMouseX = event.x
                lastMouseY = event.y
            
            ELSE IF event.type == MOUSE_UP:
                mouseDown = FALSE
            
            ELSE IF event.type == MOUSE_MOVE AND mouseDown:
                deltaX = event.x - lastMouseX
                deltaY = event.y - lastMouseY
                
                IF event.button == LEFT_BUTTON:
                    // Rotate
                    HandleMouseRotate(camera, deltaX, deltaY, 0.01)
                    renderer.needsUpdate = TRUE
                
                ELSE IF event.button == MIDDLE_BUTTON:
                    // Pan
                    HandleMousePan(camera, deltaX, deltaY, 0.01)
                    renderer.needsUpdate = TRUE
                
                lastMouseX = event.x
                lastMouseY = event.y
            
            ELSE IF event.type == MOUSE_WHEEL:
                // Zoom
                HandleMouseZoom(camera, event.delta, 0.1)
                renderer.needsUpdate = TRUE
            
            ELSE IF event.type == KEY_PRESS:
                IF event.key == SPACE:
                    // Play/pause time sequence
                    IF dataset != NULL AND dataset.numTimePoints > 1:
                        ToggleTimePlayback(dataset, renderer)
                
                ELSE IF event.key == LEFT_ARROW:
                    // Previous time point
                    IF dataset != NULL:
                        OnTimelineChanged(dataset, 
                                        dataset.currentTimeIndex - 1, 
                                        renderer)
                
                ELSE IF event.key == RIGHT_ARROW:
                    // Next time point
                    IF dataset != NULL:
                        OnTimelineChanged(dataset, 
                                        dataset.currentTimeIndex + 1, 
                                        renderer)
                
                ELSE IF event.key == R:
                    // Reset camera
                    IF dataset != NULL:
                        camera = InitializeCamera(dataset)
                        renderer.needsUpdate = TRUE
            
            ELSE IF event.type == WINDOW_RESIZE:
                // Update viewport and camera projection
                renderer.viewport.width = event.width
                renderer.viewport.height = event.height
                UpdateCameraMatrices(camera, event.width, event.height)
                renderer.needsUpdate = TRUE
        
        // Render frame if needed
        IF dataset != NULL AND renderer.needsUpdate:
            framebuffer = RenderVolume(dataset, camera, renderer)
            DisplayFramebuffer(window, framebuffer)
            renderer.needsUpdate = FALSE
        
        // Update UI
        ui.UpdateUI()
        
        // Display frame
        window.SwapBuffers()
    
    // Cleanup
    window.Close()


================================================================================
9. PERFORMANCE OPTIMIZATIONS (for future enhancements)
================================================================================

// GPU Acceleration:
// - Port ray marching to GPU shader (GLSL/HLSL)
// - Upload volume as 3D texture
// - Render to texture then display as quad
// - Massive performance improvement (100-1000x faster)

// Octree/Brick Organization:
// - Divide volume into bricks (64x64x64 voxels)
// - Build spatial index (octree)
// - Skip empty regions during ray marching
// - Reduces memory bandwidth

// Level of Detail:
// - Generate mipmaps of volume at multiple resolutions
// - Use lower resolution for distant views
// - Adaptive sampling rate based on projected size

// Progressive Rendering:
// - Render at low quality first for interactivity
// - Progressively increase quality when camera stops
// - Maintain responsive UI during interaction

// Caching:
// - Cache transfer functions as lookup textures
// - Cache camera view when time point changes
// - Pre-compute gradients for shading


================================================================================
10. TESTING AND VALIDATION
================================================================================

FUNCTION TestVolumeLoading():
    """
    Test loading various file formats
    """
    
    testFiles = [
        "test_single_channel.tif",
        "test_multi_channel.tif",
        "test_time_series.tif",
        "test_16bit.tif"
    ]
    
    FOR EACH file IN testFiles:
        dataset = LoadImageFile(file)
        ASSERT(dataset != NULL)
        ASSERT(dataset.numTimePoints > 0)
        ASSERT(dataset.GetCurrentTimePoint().channels.LENGTH > 0)


FUNCTION TestVolumeRendering():
    """
    Test rendering pipeline
    """
    
    // Create simple test volume
    volume = NEW ImageVolume()
    volume.width = 64
    volume.height = 64
    volume.depth = 64
    volume.voxelData = ALLOCATE_3D_ARRAY(64, 64, 64)
    
    // Fill with test pattern (sphere)
    center = VECTOR3(32, 32, 32)
    radius = 20
    
    FOR z = 0 TO 63:
        FOR y = 0 TO 63:
            FOR x = 0 TO 63:
                distance = SQRT((x - center.x)^2 + 
                               (y - center.y)^2 + 
                               (z - center.z)^2)
                IF distance < radius:
                    volume.SetVoxel(x, y, z, 255)
                ELSE:
                    volume.SetVoxel(x, y, z, 0)
    
    // Create minimal dataset
    dataset = CreateTestDataset(volume)
    camera = InitializeCamera(dataset)
    renderer = NEW VolumeRenderer()
    
    // Render
    framebuffer = RenderVolume(dataset, camera, renderer)
    
    // Verify output
    ASSERT(framebuffer != NULL)
    ASSERT(HasNonZeroPixels(framebuffer))


================================================================================
END OF PSEUDOCODE
================================================================================