초저온 전자현미경(Cryo-EM) 단일 입자 3차원 구조 복원 & 1,000Hz 0-Alloc C# 비정질 얼음 박막 급속 동결 및 빔 유도 드리프트 실시간 보정 서보 엔진

using System;
using System.Runtime.CompilerServices;
using System.Numerics;

namespace ZeliDesk.Stem.CryoEm
{
    /// 

    /// Cryo-EM 1,000Hz 실시간 빔 유도 이동(BIM) 추적 및 피에조 스테이지 드리프트 보정 서보
    /// GC 무할당 스택 메모리 0-Alloc 구현체
    /// 
    public ref struct CryoEmDriftCorrectionServo
    {
        // 물리 상수
        private const float ElectronWavelengthAngstrom = 0.019687f; // 300 keV 전자의 드브로이 파장
        private const float DefocusTargetNm = 1200.0f;              // 목표 디포커스 (1.2 um)
        private const float SphericalAberrationMm = 2.7f;            // 대물렌즈 구면수차 Cs (2.7 mm)
        private const float PixelSizeAngstrom = 0.85f;               // 검출기 픽셀 환산 분해능
        
        // 피에조 스테이지 서보 상태
        public float StagePosXNm;
        public float StagePosYNm;
        public float CumulativeDriftX;
        public float CumulativeDriftY;
        public float CurrentDosePerFrame; // e-/A^2
        public float TotalAccumulatedDose;
        public int ProcessedFrameCount;

        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void Initialize(float initialXNm, float initialYNm)
        {
            StagePosXNm = initialXNm;
            StagePosYNm = initialYNm;
            CumulativeDriftX = 0.0f;
            CumulativeDriftY = 0.0f;
            CurrentDosePerFrame = 1.25f; // 프레임당 1.25 e-/A^2
            TotalAccumulatedDose = 0.0f;
            ProcessedFrameCount = 0;
        }

        /// 

        /// 1ms 제어 주기: DED 서브 윈도우 위상 상관 피크로부터 나노미터 단위 드리프트 역산 및 상쇄 구동
        /// 
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public void StepFrameCorrection(
            ReadOnlySpan currentFramePatch, 
            ReadOnlySpan refFramePatch, 
            int patchSize,
            out float dXNm, 
            out float dYNm, 
            out float piezoDriveVx, 
            out float piezoDriveVy)
        {
            // 1. 위상 상관 차동 투영 피크 모멘트 탐색 (Zero-Alloc 스택 계산)
            float sumCross = 0.0f;
            float sumGradX = 0.0f;
            float sumGradY = 0.0f;

            int len = patchSize * patchSize;
            for (int i = 1; i < len - 1; i++)
            {
                float diff = currentFramePatch[i] - refFramePatch[i];
                float gradX = currentFramePatch[i + 1] - currentFramePatch[i - 1];
                float gradY = (i + patchSize < len && i - patchSize >= 0) 
                    ? (currentFramePatch[i + patchSize] - currentFramePatch[i - patchSize]) 
                    : 0.0f;

                sumCross += diff * diff;
                sumGradX += diff * gradX;
                sumGradY += diff * gradY;
            }

            // 2. 광학 유량(Optical Flow) 기반 빔 유도 서브픽셀 편위(Pixel -> Angstrom -> Nanometer)
            float denom = sumCross + 1e-6f;
            float shiftPixelX = -(sumGradX / denom);
            float shiftPixelY = -(sumGradY / denom);

            dXNm = (shiftPixelX * PixelSizeAngstrom) * 0.1f;
            dYNm = (shiftPixelY * PixelSizeAngstrom) * 0.1f;

            // 이상치 리미터 (프레임간 최대 5.0nm 제한)
            dXNm = MathF.Max(-5.0f, MathF.Min(5.0f, dXNm));
            dYNm = MathF.Max(-5.0f, MathF.Min(5.0f, dYNm));

            CumulativeDriftX += dXNm;
            CumulativeDriftY += dYNm;
            TotalAccumulatedDose += CurrentDosePerFrame;
            ProcessedFrameCount++;

            // 3. 1,000Hz 0-Alloc 능동 피에조 구동 전압 생성 (전압 감도 20 nm/V)
            const float PiezoSensitivityNmPerV = 20.0f;
            piezoDriveVx = -CumulativeDriftX / PiezoSensitivityNmPerV;
            piezoDriveVy = -CumulativeDriftY / PiezoSensitivityNmPerV;

            // 스테이지 피드백 갱신
            StagePosXNm += dXNm;
            StagePosYNm += dYNm;
        }

        /// 

        /// 공간주파수 k(1/A)에 대한 CTF 위상값 및 손상 감쇠 가중치 계산
        /// 
        [MethodImpl(MethodImplOptions.AggressiveInlining)]
        public readonly float ComputeExposureWeight(float spatialFreqInvA)
        {
            // 임계 도스 모델 Ne(k) = 0.245 * k^(-1.665) + 2.81
            float k = MathF.Max(0.01f, spatialFreqInvA);
            float criticalDose = 0.245f * MathF.Pow(k, -1.665f) + 2.81f;
            return MathF.Exp(-TotalAccumulatedDose / (2.0f * criticalDose));
        }
    }
}