# Copyright (c) DataLab Platform Developers, BSD 3-Clause license, see LICENSE file.""".. Image Processing Algorithms (see parent package :mod:`cdl.algorithms`)"""# pylint: disable=invalid-name # Allows short reference names like x, y, ...from__future__importannotationsfromtypingimportLiteralimportnumpyasnpimportscipy.ndimageasspiimportscipy.spatialassptfromnumpyimportmafromskimageimportexposure,feature,measure,transform# MARK: Level adjustment ---------------------------------------------------------------
[docs]defscale_data_to_min_max(data:np.ndarray,zmin:float|int,zmax:float|int)->np.ndarray:"""Scale array `data` to fit [zmin, zmax] dynamic range Args: data: Input data zmin: Minimum value of output data zmax: Maximum value of output data Returns: Scaled data """dmin=data.min()dmax=data.max()fdata=np.array(data,dtype=float)fdata-=dminfdata*=float(zmax-zmin)/(dmax-dmin)fdata+=float(zmin)returnnp.array(fdata,data.dtype)
[docs]defnormalize(data:np.ndarray,parameter:Literal["maximum","amplitude","area","energy","rms"]="maximum",)->np.ndarray:"""Normalize input array to a given parameter. Args: data: Input data parameter: Normalization parameter (default: "maximum") Returns: Normalized array """ifparameter=="maximum":returnscale_data_to_min_max(data,data.min()/data.max(),1.0)ifparameter=="amplitude":returnscale_data_to_min_max(data,0.0,1.0)fdata=np.array(data,dtype=float)ifparameter=="area":returnfdata/fdata.sum()ifparameter=="energy":returnfdata/np.sqrt(np.sum(fdata*fdata.conjugate()))ifparameter=="rms":returnfdata/np.sqrt(np.mean(fdata*fdata.conjugate()))raiseValueError(f"Unsupported parameter {parameter}")
[docs]deffft2d(z:np.ndarray,shift:bool=True)->np.ndarray:"""Compute FFT of complex array `z` Args: z: Input data shift: Shift zero frequency to center (default: True) Returns: FFT of input data """z1=np.fft.fft2(z)ifshift:z1=np.fft.fftshift(z1)returnz1
[docs]defifft2d(z:np.ndarray,shift:bool=True)->np.ndarray:"""Compute inverse FFT of complex array `z` Args: z: Input data shift: Shift zero frequency to center (default: True) Returns: Inverse FFT of input data """ifshift:z=np.fft.ifftshift(z)z1=np.fft.ifft2(z)returnz1
[docs]defmagnitude_spectrum(z:np.ndarray,log_scale:bool=False)->np.ndarray:"""Compute magnitude spectrum of complex array `z` Args: z: Input data log_scale: Use log scale (default: False) Returns: Magnitude spectrum of input data """z1=np.abs(fft2d(z))iflog_scale:z1=20*np.log10(z1.clip(1e-10))returnz1
[docs]defphase_spectrum(z:np.ndarray)->np.ndarray:"""Compute phase spectrum of complex array `z` Args: z: Input data Returns: Phase spectrum of input data (in degrees) """returnnp.rad2deg(np.angle(fft2d(z)))
[docs]defpsd(z:np.ndarray,log_scale:bool=False)->np.ndarray:"""Compute power spectral density of complex array `z` Args: z: Input data log_scale: Use log scale (default: False) Returns: Power spectral density of input data """z1=np.abs(fft2d(z))**2iflog_scale:z1=10*np.log10(z1.clip(1e-10))returnz1
[docs]defbinning(data:np.ndarray,sx:int,sy:int,operation:Literal["sum","average","median","min","max"],dtype=None,)->np.ndarray:"""Perform image pixel binning Args: data: Input data sx: Binning size along x (number of pixels to bin together) sy: Binning size along y (number of pixels to bin together) operation: Binning operation dtype: Output data type (default: None, i.e. same as input) Returns: Binned data """ny,nx=data.shapeshape=(ny//sy,sy,nx//sx,sx)try:bdata=data[:ny-ny%sy,:nx-nx%sx].reshape(shape)exceptValueErroraserr:raiseValueError("Binning is not a multiple of image dimensions")fromerrifoperation=="sum":bdata=np.array(bdata,dtype=float).sum(axis=(-1,1))elifoperation=="average":bdata=bdata.mean(axis=(-1,1))elifoperation=="median":bdata=ma.median(bdata,axis=(-1,1))elifoperation=="min":bdata=bdata.min(axis=(-1,1))elifoperation=="max":bdata=bdata.max(axis=(-1,1))else:valid=", ".join(BINNING_OPERATIONS)raiseValueError(f"Invalid operation {operation} (valid values: {valid})")returnnp.array(bdata,dtype=data.dtypeifdtypeisNoneelsenp.dtype(dtype))
[docs]defflatfield(rawdata:np.ndarray,flatdata:np.ndarray,threshold:float|None=None)->np.ndarray:"""Compute flat-field correction Args: rawdata: Raw data flatdata: Flat-field data threshold: Threshold for flat-field correction (default: None) Returns: Flat-field corrected data """dtemp=np.array(rawdata,dtype=float,copy=True)*flatdata.mean()dunif=np.array(flatdata,dtype=float,copy=True)dunif[dunif==0]=1.0dcorr_all=np.array(dtemp/dunif,dtype=rawdata.dtype)dcorr=np.array(rawdata,copy=True)dcorr[rawdata>threshold]=dcorr_all[rawdata>threshold]returndcorr
[docs]defget_centroid_fourier(data:np.ndarray)->tuple[float,float]:"""Return image centroid using Fourier algorithm Args: data: Input data Returns: Centroid coordinates (row, col) """# Fourier transform method as discussed by Weisshaar et al.# (http://www.mnd-umwelttechnik.fh-wiesbaden.de/pig/weisshaar_u5.pdf)rows,cols=data.shapeifrows==1orcols==1:return0,0i=np.arange(0,rows).reshape(1,rows)sin_a=np.sin((i-1)*2*np.pi/(rows-1)).Tcos_a=np.cos((i-1)*2*np.pi/(rows-1)).Tj=np.arange(0,cols).reshape(cols,1)sin_b=np.sin((j-1)*2*np.pi/(cols-1)).Tcos_b=np.cos((j-1)*2*np.pi/(cols-1)).Ta=(cos_a*data).sum()b=(sin_a*data).sum()c=(data*cos_b).sum()d=(data*sin_b).sum()rphi=(0ifb>0else2*np.pi)ifa>0elsenp.picphi=(0ifd>0else2*np.pi)ifc>0elsenp.piifa*c==0.0:return0,0row=(np.arctan(b/a)+rphi)*(rows-1)/(2*np.pi)+1col=(np.arctan(d/c)+cphi)*(cols-1)/(2*np.pi)+1row=np.nanifrowisnp.ma.maskedelserowcol=np.nanifcolisnp.ma.maskedelsecolreturnrow,col
[docs]defget_absolute_level(data:np.ndarray,level:float)->float:"""Return absolute level Args: data: Input data level: Relative level (0.0 to 1.0) Returns: Absolute level """ifnotisinstance(level,float)orlevel<0.0orlevel>1.0:raiseValueError("Level must be a float between 0. and 1.")return(float(np.nanmin(data))+float(np.nanmax(data)))*level
[docs]defget_enclosing_circle(data:np.ndarray,level:float=0.5)->tuple[int,int,float]:"""Return (x, y, radius) for the circle contour enclosing image values above threshold relative level (.5 means FWHM) Args: data: Input data level: Relative level (default: 0.5) Returns: A tuple (x, y, radius) Raises: ValueError: No contour was found """data_th=data.copy()data_th[data<=get_absolute_level(data,level)]=0.0contours=measure.find_contours(data_th)model=measure.CircleModel()result=Nonemax_radius=1.0forcontourincontours:ifmodel.estimate(contour):yc,xc,radius=model.paramsifradius>max_radius:result=(int(xc),int(yc),radius)max_radius=radiusifresultisNone:raiseValueError("No contour was found")returnresult
[docs]defget_radial_profile(data:np.ndarray,center:tuple[int,int])->tuple[np.ndarray,np.ndarray]:"""Return radial profile of image data Args: data: Input data (2D array) center: Coordinates of the center of the profile (x, y) Returns: Radial profile (X, Y) where X is the distance from the center (1D array) and Y is the average value of pixels at this distance (1D array) """y,x=np.indices((data.shape))# Get the indices of pixelsx0,y0=centerr=np.sqrt((x-x0)**2+(y-y0)**2)# Calculate distance to the centerr=r.astype(int)# Average over the same distancetbin=np.bincount(r.ravel(),data.ravel())# Sum of pixel values at each distancenr=np.bincount(r.ravel())# Number of pixels at each distanceyprofile=tbin/nr# this is the half radial profile# Let's mirror it to get the full radial profile (the first element is the center)yprofile=np.concatenate((yprofile[::-1],yprofile[1:]))# The x axis is the distance from the center (0 is the center)xprofile=np.arange(len(yprofile))-len(yprofile)//2returnxprofile,yprofile
[docs]defdistance_matrix(coords:list)->np.ndarray:"""Return distance matrix from coords Args: coords: List of coordinates Returns: Distance matrix """returnnp.triu(spt.distance.cdist(coords,coords,"euclidean"))
[docs]defget_2d_peaks_coords(data:np.ndarray,size:int|None=None,level:float=0.5)->np.ndarray:"""Detect peaks in image data, return coordinates. If neighborhoods size is None, default value is the highest value between 50 pixels and the 1/40th of the smallest image dimension. Detection threshold level is relative to difference between data maximum and minimum values. Args: data: Input data size: Neighborhood size (default: None) level: Relative level (default: 0.5) Returns: Coordinates of peaks """ifsizeisNone:size=max(min(data.shape)//40,50)data_max=spi.maximum_filter(data,size)data_min=spi.minimum_filter(data,size)data_diff=data_max-data_mindiff=(data_max-data_min)>get_absolute_level(data_diff,level)maxima=data==data_maxmaxima[diff==0]=0labeled,_num_objects=spi.label(maxima)slices=spi.find_objects(labeled)coords=[]fordy,dxinslices:x_center=int(0.5*(dx.start+dx.stop-1))y_center=int(0.5*(dy.start+dy.stop-1))coords.append((x_center,y_center))iflen(coords)>1:# Eventually removing duplicatesdist=distance_matrix(coords)forindexinreversed(np.unique(np.where((dist<size)&(dist>0))[1])):coords.pop(index)returnnp.array(coords)
[docs]defget_contour_shapes(data:np.ndarray|ma.MaskedArray,shape:Literal["circle","ellipse","polygon"]="ellipse",level:float=0.5,)->np.ndarray:"""Find iso-valued contours in a 2D array, above relative level (.5 means FWHM), then fit contours with shape ('ellipse' or 'circle') Args: data: Input data shape: Shape to fit. Default is 'ellipse' level: Relative level (default: 0.5) Returns: Coordinates of shapes """# pylint: disable=too-many-localsassertshapein("circle","ellipse","polygon")contours=measure.find_contours(data,level=get_absolute_level(data,level))coords=[]forcontourincontours:# `contour` is a (N, 2) array (rows, cols): we need to check if all those# coordinates are masked: if so, we skip this contourifisinstance(data,ma.MaskedArray)andnp.all(data.mask[contour[:,0].astype(int),contour[:,1].astype(int)]):continueifshape=="circle":model=measure.CircleModel()ifmodel.estimate(contour):yc,xc,r=model.paramsifr<=1.0:continuecoords.append([xc,yc,r])elifshape=="ellipse":model=measure.EllipseModel()ifmodel.estimate(contour):yc,xc,b,a,theta=model.paramsifa<=1.0orb<=1.0:continuecoords.append([xc,yc,a,b,theta])elifshape=="polygon":# `contour` is a (N, 2) array (rows, cols): we need to convert it# to a list of x, y coordinates flattened in a single listcoords.append(contour[:,::-1].flatten())else:raiseNotImplementedError(f"Invalid contour model {model}")ifshape=="polygon":# `coords` is a list of arrays of shape (N, 2) where N is the number of points# that can vary from one array to another, so we need to padd with NaNs each# array to get a regular array:max_len=max(coord.shape[0]forcoordincoords)arr=np.full((len(coords),max_len),np.nan)fori_row,coordinenumerate(coords):arr[i_row,:coord.shape[0]]=coordreturnarrreturnnp.array(coords)
[docs]defget_hough_circle_peaks(data:np.ndarray,min_radius:float|None=None,max_radius:float|None=None,nb_radius:int|None=None,min_distance:int=1,)->np.ndarray:"""Detect peaks in image from circle Hough transform, return circle coordinates. Args: data: Input data min_radius: Minimum radius (default: None) max_radius: Maximum radius (default: None) nb_radius: Number of radii (default: None) min_distance: Minimum distance between circles (default: 1) Returns: Coordinates of circles """assertmin_radiusisnotNoneandmax_radiusisnotNoneandmax_radius>min_radiusifnb_radiusisNone:nb_radius=max_radius-min_radius+1hough_radii=np.arange(min_radius,max_radius+1,(max_radius-min_radius+1)//nb_radius)hough_res=transform.hough_circle(data,hough_radii)_accums,cx,cy,radii=transform.hough_circle_peaks(hough_res,hough_radii,min_xdistance=min_distance,min_ydistance=min_distance)returnnp.vstack([cx,cy,radii]).T
[docs]deffind_blobs_dog(data:np.ndarray,min_sigma:float=1,max_sigma:float=30,overlap:float=0.5,threshold_rel:float=0.2,exclude_border:bool=True,)->np.ndarray:""" Finds blobs in the given grayscale image using the Difference of Gaussians (DoG) method. Args: data: The grayscale input image. min_sigma: The minimum blob radius in pixels. max_sigma: The maximum blob radius in pixels. overlap: The minimum overlap between two blobs in pixels. For instance, if two blobs are detected with radii of 10 and 12 respectively, and the ``overlap`` is set to 0.5, then the area of the smaller blob will be ignored and only the area of the larger blob will be returned. threshold_rel: The absolute lower bound for scale space maxima. Local maxima smaller than ``threshold_rel`` are ignored. Reduce this to detect blobs with less intensities. exclude_border: If ``True``, exclude blobs from detection if they are too close to the border of the image. Border size is ``min_sigma``. Returns: Coordinates of blobs """# Use scikit-image's Difference of Gaussians (DoG) methodblobs=feature.blob_dog(data,min_sigma=min_sigma,max_sigma=max_sigma,overlap=overlap,threshold_rel=threshold_rel,exclude_border=exclude_border,)return__blobs_to_coords(blobs)
[docs]deffind_blobs_doh(data:np.ndarray,min_sigma:float=1,max_sigma:float=30,overlap:float=0.5,log_scale:bool=False,threshold_rel:float=0.2,)->np.ndarray:""" Finds blobs in the given grayscale image using the Determinant of Hessian (DoH) method. Args: data: The grayscale input image. min_sigma: The minimum blob radius in pixels. max_sigma: The maximum blob radius in pixels. overlap: The minimum overlap between two blobs in pixels. For instance, if two blobs are detected with radii of 10 and 12 respectively, and the ``overlap`` is set to 0.5, then the area of the smaller blob will be ignored and only the area of the larger blob will be returned. log_scale: If ``True``, the radius of each blob is returned as ``sqrt(sigma)`` for each detected blob. threshold_rel: The absolute lower bound for scale space maxima. Local maxima smaller than ``threshold_rel`` are ignored. Reduce this to detect blobs with less intensities. Returns: Coordinates of blobs """# Use scikit-image's Determinant of Hessian (DoH) method to detect blobsblobs=feature.blob_doh(data,min_sigma=min_sigma,max_sigma=max_sigma,num_sigma=int(max_sigma-min_sigma+1),threshold=None,threshold_rel=threshold_rel,overlap=overlap,log_scale=log_scale,)return__blobs_to_coords(blobs)
[docs]deffind_blobs_log(data:np.ndarray,min_sigma:float=1,max_sigma:float=30,overlap:float=0.5,log_scale:bool=False,threshold_rel:float=0.2,exclude_border:bool=True,)->np.ndarray:"""Finds blobs in the given grayscale image using the Laplacian of Gaussian (LoG) method. Args: data: The grayscale input image. min_sigma: The minimum blob radius in pixels. max_sigma: The maximum blob radius in pixels. overlap: The minimum overlap between two blobs in pixels. For instance, if two blobs are detected with radii of 10 and 12 respectively, and the ``overlap`` is set to 0.5, then the area of the smaller blob will be ignored and only the area of the larger blob will be returned. log_scale: If ``True``, the radius of each blob is returned as ``sqrt(sigma)`` for each detected blob. threshold_rel: The absolute lower bound for scale space maxima. Local maxima smaller than ``threshold_rel`` are ignored. Reduce this to detect blobs with less intensities. exclude_border: If ``True``, exclude blobs from detection if they are too close to the border of the image. Border size is ``min_sigma``. Returns: Coordinates of blobs """# Use scikit-image's Laplacian of Gaussian (LoG) method to detect blobsblobs=feature.blob_log(data,min_sigma=min_sigma,max_sigma=max_sigma,num_sigma=int(max_sigma-min_sigma+1),threshold=None,threshold_rel=threshold_rel,overlap=overlap,log_scale=log_scale,exclude_border=exclude_border,)return__blobs_to_coords(blobs)
[docs]defremove_overlapping_disks(coords:np.ndarray)->np.ndarray:"""Remove overlapping disks among coordinates Args: coords: The coordinates of the disks Returns: The coordinates of the disks with overlapping disks removed """# Get the radii of each disk from the coordinatesradii=coords[:,2]# Calculate the distance between the center of each pair of disksdist=np.sqrt(np.sum((coords[:,None,:2]-coords[:,:2])**2,axis=-1))# Create a boolean mask where the distance between the centers# is less than the sum of the radiimask=dist<(radii[:,None]+radii)# Find the indices of overlapping disksoverlapping_indices=np.argwhere(mask)# Remove the smaller disk from each overlapping pairfori,jinoverlapping_indices:ifi!=j:ifradii[i]<radii[j]:coords[i]=[np.nan,np.nan,np.nan]else:coords[j]=[np.nan,np.nan,np.nan]# Remove rows with NaN valuescoords=coords[~np.isnan(coords).any(axis=1)]returncoords
[docs]deffind_blobs_opencv(data:np.ndarray,min_threshold:float|None=None,max_threshold:float|None=None,min_repeatability:int|None=None,min_dist_between_blobs:float|None=None,filter_by_color:bool|None=None,blob_color:int|None=None,filter_by_area:bool|None=None,min_area:float|None=None,max_area:float|None=None,filter_by_circularity:bool|None=None,min_circularity:float|None=None,max_circularity:float|None=None,filter_by_inertia:bool|None=None,min_inertia_ratio:float|None=None,max_inertia_ratio:float|None=None,filter_by_convexity:bool|None=None,min_convexity:float|None=None,max_convexity:float|None=None,)->np.ndarray:""" Finds blobs in the given grayscale image using OpenCV's SimpleBlobDetector. Args: data: The grayscale input image. min_threshold: The minimum blob intensity. max_threshold: The maximum blob intensity. min_repeatability: The minimum number of times a blob is detected before it is reported. min_dist_between_blobs: The minimum distance between blobs. filter_by_color: If ``True``, blobs are filtered by color. blob_color: The color of the blobs to filter by. filter_by_area: If ``True``, blobs are filtered by area. min_area: The minimum blob area. max_area: The maximum blob area. filter_by_circularity: If ``True``, blobs are filtered by circularity. min_circularity: The minimum blob circularity. max_circularity: The maximum blob circularity. filter_by_inertia: If ``True``, blobs are filtered by inertia. min_inertia_ratio: The minimum blob inertia ratio. max_inertia_ratio: The maximum blob inertia ratio. filter_by_convexity: If ``True``, blobs are filtered by convexity. min_convexity: The minimum blob convexity. max_convexity: The maximum blob convexity. Returns: Coordinates of blobs """# Note:# Importing OpenCV inside the function in order to eventually raise an ImportError# when the function is called and OpenCV is not installed. This error will be# handled by DataLab and the user will be informed that OpenCV is required to use# this function.importcv2# pylint: disable=import-outside-toplevelparams=cv2.SimpleBlobDetector_Params()ifmin_thresholdisnotNone:params.minThreshold=min_thresholdifmax_thresholdisnotNone:params.maxThreshold=max_thresholdifmin_repeatabilityisnotNone:params.minRepeatability=min_repeatabilityifmin_dist_between_blobsisnotNone:params.minDistBetweenBlobs=min_dist_between_blobsiffilter_by_colorisnotNone:params.filterByColor=filter_by_colorifblob_colorisnotNone:params.blobColor=blob_coloriffilter_by_areaisnotNone:params.filterByArea=filter_by_areaifmin_areaisnotNone:params.minArea=min_areaifmax_areaisnotNone:params.maxArea=max_areaiffilter_by_circularityisnotNone:params.filterByCircularity=filter_by_circularityifmin_circularityisnotNone:params.minCircularity=min_circularityifmax_circularityisnotNone:params.maxCircularity=max_circularityiffilter_by_inertiaisnotNone:params.filterByInertia=filter_by_inertiaifmin_inertia_ratioisnotNone:params.minInertiaRatio=min_inertia_ratioifmax_inertia_ratioisnotNone:params.maxInertiaRatio=max_inertia_ratioiffilter_by_convexityisnotNone:params.filterByConvexity=filter_by_convexityifmin_convexityisnotNone:params.minConvexity=min_convexityifmax_convexityisnotNone:params.maxConvexity=max_convexitydetector=cv2.SimpleBlobDetector_create(params)image=exposure.rescale_intensity(data,out_range=np.uint8)keypoints=detector.detect(image)ifkeypoints:coords=cv2.KeyPoint_convert(keypoints)radii=0.5*np.array([kp.sizeforkpinkeypoints])blobs=np.vstack([coords[:,1],coords[:,0],radii]).Tblobs=remove_overlapping_disks(blobs)else:blobs=np.array([]).reshape((0,3))return__blobs_to_coords(blobs)