"""Python adaptation of OpenBreast algorithms (Said Pertuz / F. Torres, BSD-2).
See third_party/OpenBreast/LICENSE and PORT_NOTES.md. Not a validated MATLAB clone.
Working coordinates: breast/chest wall on LEFT; arrays in [row, column] order.
"""
import numpy as np
from scipy import ndimage
from scipy.optimize import curve_fit
from scipy.signal import find_peaks
from scipy.signal.windows import gaussian
from skimage import feature, measure, morphology, filters

VERSION='OpenBreast-Python-adaptation-1.0'

def largest(mask):
    labels,n=ndimage.label(mask,structure=np.ones((3,3)))
    if not n: raise ValueError('유방 영역을 찾지 못했습니다.')
    counts=np.bincount(labels.ravel()); counts[0]=0
    return labels==counts.argmax()

def breast_on_right(im):
    tmp=im.copy(); tmp[tmp>.95*im.max()]=0
    if tmp.sum()==0:tmp=im
    mid=im.shape[1]//2
    return bool(tmp[:,mid:].sum()>tmp[:,:mid].sum())

def foreground_threshold(im):
    """OpenBreast histogram Gaussian fit; named Otsu fallback if fit fails."""
    if not np.isfinite(im).all() or np.ptp(im)<1e-6: raise ValueError('영상 대비가 없거나 비정상 pixel이 있습니다.')
    x=im.ravel(); xi=np.linspace(float(x.min()),float(x.max()),1000)
    hist,_=np.histogram(x,bins=np.r_[xi,xi[-1]+1e-6])
    n=ndimage.convolve1d(hist.astype(float),gaussian(25,4.8),mode='constant')
    n/=max(n.max(),1)
    xsup=min(float(x.mean()),max(float(np.percentile(x,30)),.2))
    peaks,_=find_peaks(n*(xi<=xsup),height=.35)
    method='OpenBreast histogram fit'
    try:
        if len(peaks):
            select=(n>.35)&(xi<xsup)
            if len(peaks)>1:
                start=peaks.min()+np.argmin(n[peaks.min():peaks.max()+1]); select &= xi>=xi[start]
            if select.sum()<4: raise ValueError('Insufficient histogram fit points')
            def gauss(t,a,b,c): return a*np.exp(-((t-b)/c)**2)
            t=xi[select]; y=n[select]
            params,_=curve_fit(gauss,t,y,p0=[y.max(),t[y.argmax()],max(np.std(t),.005)],
                              bounds=([.05,float(x.min()),1e-5],[3,float(x.max()),1]),maxfev=4000)
            a,b,c=params; threshold=float(b+np.sqrt(c*c*np.log(a/.05)))
        else:
            threshold=float(xi[np.flatnonzero(n<.05*n.max())[0]])
        if not x.min()<threshold<x.max(): raise ValueError('Invalid histogram threshold')
    except (ValueError,RuntimeError,FloatingPointError):
        threshold=float(filters.threshold_otsu(im)); method='Otsu fallback (histogram fit failure)'
    return threshold,method


def foreground(im):
    threshold,method=foreground_threshold(im)
    mask=im>=max(threshold,0)
    mask[0,:]=False; mask[-1,:]=False
    cleaned=ndimage.binary_erosion(mask,structure=np.ones((5,5)))
    if not cleaned.any(): raise ValueError('윤곽 추출 실패; threshold 또는 mask를 수정하세요.')
    cleaned=largest(cleaned)
    cleaned=ndimage.binary_dilation(cleaned,structure=np.ones((5,5)))
    cleaned=ndimage.binary_closing(cleaned,structure=np.ones((5,5)))
    cleaned=ndimage.binary_fill_holes(cleaned)
    fraction=cleaned.mean()
    warnings=[]
    if fraction<.03 or fraction>.95: warnings.append('유방 면적 비율이 비정상입니다. 윤곽을 확인하세요.')
    return cleaned,threshold,method,warnings

def pectoral_line(im,mask,processor=None):
    """Adapt OpenBreast segChest Canny + restricted rho/theta voting.
    Gate output to avoid silently deleting breast when evidence is weak.
    """
    h,w=im.shape; ys,xs=np.where(mask)
    ymax=max(5,min(h,int(.6*(ys.max()+1))))
    row_extents=[]
    for y in range(max(1,ys.min()),ymax):
        xx=np.flatnonzero(mask[y])
        if len(xx): row_extents.append(xx.max()+1)
    xmax=max(5,min(w,int(np.percentile(row_extents,20)))) if row_extents else max(5,w//2)
    smooth=processor.pectoral_smooth(im) if processor else ndimage.gaussian_filter(morphology.dilation(im,morphology.disk(8)),1)
    edges=feature.canny(smooth,sigma=2)[:ymax,:xmax]
    yy,xx=np.indices(edges.shape); edges[yy>ymax-(ymax-1)*xx/max(xmax-1,1)]=False
    y,x=np.where(edges)
    if len(x)<20: return None,'대흉근 선 근거 부족: 수동 확인 필요'
    theta=np.linspace(np.deg2rad(20),np.deg2rad(45),128)
    rho=np.linspace(1,min(edges.shape),128)
    accum=np.zeros((128,128),int)
    for j,t in enumerate(theta):
        r=x*np.cos(t)+y*np.sin(t)
        accum[:,j]=np.histogram(r,bins=np.r_[rho,rho[-1]+(rho[1]-rho[0])])[0]
    i,j=np.unravel_index(accum.argmax(),accum.shape)
    if accum[i,j]<max(12,.035*len(x)): return None,'대흉근 선 신뢰도 부족: 수동 확인 필요'
    slope=-np.cos(theta[j])/np.sin(theta[j]); intercept=rho[i]/np.sin(theta[j])
    yy,xx=np.indices(mask.shape)
    removed=mask & (yy<intercept+slope*xx)
    fraction=removed.sum()/mask.sum()
    if fraction<.002 or fraction>.35: return None,'대흉근 제외 면적이 비정상: 수동 확인 필요'
    return {'slope':float(slope),'intercept_px':float(intercept),'removed_fraction':float(fraction)},None

def apply_line(mask,line):
    if line is None: return mask.copy()
    yy,xx=np.indices(mask.shape)
    result=mask & (yy>=line['intercept_px']+line['slope']*xx)
    if not result.any(): raise ValueError('대흉근 선이 유방 전체를 제거했습니다.')
    return result

def skin_source(mask,skin_mask=None):
    """Mask whose boundary is skin only: before pectoral removal, extended to the image borders.
    foreground() leaves the mask a few pixels short of the image border, which would otherwise be
    eroded like skin. Working coordinates have the chest wall on the LEFT.
    """
    source=np.array(mask if skin_mask is None else skin_mask,dtype=bool)
    # Everything between the chest-wall image edge and the breast in the same row is not skin.
    source=np.logical_or.accumulate(source[:,::-1],axis=1)[:,::-1]
    margin=min(3,source.shape[0]//2)
    if margin: source[:margin]=source[margin]; source[-margin:]=source[-margin-1]
    return source

def dense_threshold(mag_threshold,lo,hi,offset=0.):
    """Automatic MAG threshold plus a manual offset, kept inside the breast intensity range.
    Returns (automatic, used)."""
    auto=float(mag_threshold)
    return auto,float(min(max(auto+offset,lo),hi))

def mag_density(im,mask,step_mm,skin_gap_mm=4.,min_area_mm2=16.,skin_mask=None,offset=0.):
    """Port mseg/pprocess: 25 levels, MAG Gaussian smoothing, erode and remove islands.
    Fraction returned as percent, NOT BI-RADS or volumetric density.
    skin_mask: breast mask before pectoral removal; the skin gap is eroded from it so that
    neither the chest-wall image edge nor the pectoral line gets an extra margin.
    offset: see dense_threshold; the default reproduces the original MAG threshold.
    """
    mask=np.asarray(mask,bool)
    if not mask.any(): raise ValueError('빈 유방 mask')
    lo,hi=float(im[mask].min()),float(im[mask].max())
    if hi-lo<1e-6: raise ValueError('유방 내부 대비가 없어 밀도를 계산할 수 없습니다.')
    levels=np.linspace(lo,hi,25)
    area=np.array([np.count_nonzero((im>=v)&mask) for v in levels],float)
    # MATLAB gausswin(5), alpha=2.5 => standard deviation 0.8 pixels.
    gradient=np.convolve(np.diff(area),gaussian(5,.8),'same')
    index=int(np.argmin(gradient))
    auto,threshold=dense_threshold(levels[index+1],lo,hi,offset)
    dense=(im>=threshold)&mask
    warnings=[]
    if step_mm is None:
        raise ValueError('MAG 밀도 계산에는 mm pixel spacing이 필요합니다. 태그 또는 수동 값을 입력하세요.')
    diameter=max(1,int(np.floor(2*skin_gap_mm/step_mm+.5)))
    # min-filter is equivalent to erosion by a square footprint and avoids huge footprints.
    # cval=1: outside the image counts as breast (as MATLAB imerode), so image edges are not eroded.
    eroded=(ndimage.minimum_filter(skin_source(mask,skin_mask).astype('uint8'),size=diameter,mode='constant',cval=1)>0)&mask
    if not eroded.any(): raise ValueError('피부 제외 후 분석 영역이 없습니다. 윤곽 / spacing / skin gap 확인')
    minimum=max(1,int(np.floor(min_area_mm2/(step_mm**2)+.5)))
    dense=morphology.remove_small_objects(dense&eroded,min_size=minimum,connectivity=2)
    pd=float(100*dense.sum()/mask.sum())
    if pd<.1 or pd>90: warnings.append('극단적인 MAG 밀도 값: overlay와 전처리를 검토하세요.')
    return dense,pd,{'threshold':threshold,'automatic_threshold':auto,'threshold_offset':offset,
                      'skin_gap_mm':skin_gap_mm,'minimum_component_mm2':min_area_mm2,'method':VERSION,'denominator':'whole breast mask, including skin gap',
                      'levels':25,'area_curve':area.tolist(),'warnings':warnings}

def analyze(im,ismlo,step_mm,orientation='auto',manual_threshold=None,manual_line=None,skip_pectoral=False,skin_gap_mm=4.,inferior_cut_fraction=None,processor=None,dense_offset=.07):
    flip=breast_on_right(im) if orientation=='auto' else orientation=='right'
    canonical=np.fliplr(im) if flip else im
    if manual_threshold is not None:
        mask=ndimage.binary_fill_holes(largest(canonical>=manual_threshold));threshold=manual_threshold;method='manual threshold';warnings=[]
    else:
        mask,threshold,method,warnings=processor.foreground(canonical) if processor else foreground(canonical)
    if inferior_cut_fraction is not None:
        cutoff=max(1,min(mask.shape[0],int(round(mask.shape[0]*inferior_cut_fraction))))
        mask[cutoff:,:]=False
        if not mask.any():raise ValueError('하부 절단 후 유방 영역이 없습니다.')
    elif ismlo and mask[-max(4,int(.03*mask.shape[0])):,:].any():
        warnings.append('MLO 윤곽이 영상 하단에 이어집니다. 하부 피부 주름 포함 여부를 확인하고 필요 시 하부 절단을 설정하세요.')
    line=None;skin_mask=mask
    if ismlo and not skip_pectoral:
        if manual_line is not None: line=manual_line
        else:
            line,message=pectoral_line(canonical,mask,processor)
            if message: warnings.append(message)
        mask=apply_line(mask,line)
    if ismlo and skip_pectoral: warnings.append('MLO 대흉근 제외 생략: 밀도·면적에 영향 가능')
    try:
        dense,pd,detail=(processor.mag_density if processor else mag_density)(canonical,mask,step_mm,skin_gap_mm,skin_mask=skin_mask,offset=dense_offset)
        warnings+=detail.pop('warnings')
    except ValueError as exc:
        dense=np.zeros_like(mask);pd=None;detail={'method':VERSION,'error':str(exc)};warnings.append(str(exc))
    return {'breast_mask':np.fliplr(mask) if flip else mask,'dense_mask':np.fliplr(dense) if flip else dense,
            'density_percent':pd,'threshold':threshold,'foreground_method':method,'pectoral_line':line,
            'flipped_for_analysis':flip,'mag':detail,'warnings':warnings}