Mercurial > repos > perssond > unmicst
diff toolbox/imtools.py @ 1:74fe58ff55a5 draft default tip
planemo upload for repository https://github.com/HMS-IDAC/UnMicst commit e14f76a8803cab0013c6dbe809bc81d7667f2ab9
author | goeckslab |
---|---|
date | Wed, 07 Sep 2022 23:10:14 +0000 |
parents | 6bec4fef6b2e |
children |
line wrap: on
line diff
--- a/toolbox/imtools.py Fri Mar 12 00:17:29 2021 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,312 +0,0 @@ -import matplotlib.pyplot as plt -import tifffile -import os -import numpy as np -from skimage import io as skio -from scipy.ndimage import * -from skimage.morphology import * -from skimage.transform import resize - -def tifread(path): - return tifffile.imread(path) - -def tifwrite(I,path): - tifffile.imsave(path, I) - -def imshow(I,**kwargs): - if not kwargs: - plt.imshow(I,cmap='gray') - else: - plt.imshow(I,**kwargs) - - plt.axis('off') - plt.show() - -def imshowlist(L,**kwargs): - n = len(L) - for i in range(n): - plt.subplot(1, n, i+1) - if not kwargs: - plt.imshow(L[i],cmap='gray') - else: - plt.imshow(L[i],**kwargs) - plt.axis('off') - plt.show() - -def imread(path): - return skio.imread(path) - -def imwrite(I,path): - skio.imsave(path,I) - -def im2double(I): - if I.dtype == 'uint16': - return I.astype('float64')/65535 - elif I.dtype == 'uint8': - return I.astype('float64')/255 - elif I.dtype == 'float32': - return I.astype('float64') - elif I.dtype == 'float64': - return I - else: - print('returned original image type: ', I.dtype) - return I - -def size(I): - return list(I.shape) - -def imresizeDouble(I,sizeOut): # input and output are double - return resize(I,(sizeOut[0],sizeOut[1]),mode='reflect') - -def imresize3Double(I,sizeOut): # input and output are double - return resize(I,(sizeOut[0],sizeOut[1],sizeOut[2]),mode='reflect') - -def imresizeUInt8(I,sizeOut): # input and output are UInt8 - return np.uint8(resize(I.astype(float),(sizeOut[0],sizeOut[1]),mode='reflect',order=0)) - -def imresize3UInt8(I,sizeOut): # input and output are UInt8 - return np.uint8(resize(I.astype(float),(sizeOut[0],sizeOut[1],sizeOut[2]),mode='reflect',order=0)) - -def normalize(I): - m = np.min(I) - M = np.max(I) - if M > m: - return (I-m)/(M-m) - else: - return I - -def snormalize(I): - m = np.mean(I) - s = np.std(I) - if s > 0: - return (I-m)/s - else: - return I - -def cat(a,I,J): - return np.concatenate((I,J),axis=a) - -def imerode(I,r): - return binary_erosion(I, disk(r)) - -def imdilate(I,r): - return binary_dilation(I, disk(r)) - -def imerode3(I,r): - return morphology.binary_erosion(I, ball(r)) - -def imdilate3(I,r): - return morphology.binary_dilation(I, ball(r)) - -def sphericalStructuralElement(imShape,fRadius): - if len(imShape) == 2: - return disk(fRadius,dtype=float) - if len(imShape) == 3: - return ball(fRadius,dtype=float) - -def medfilt(I,filterRadius): - return median_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius)) - -def maxfilt(I,filterRadius): - return maximum_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius)) - -def minfilt(I,filterRadius): - return minimum_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius)) - -def ptlfilt(I,percentile,filterRadius): - return percentile_filter(I,percentile,footprint=sphericalStructuralElement(I.shape,filterRadius)) - -def imgaussfilt(I,sigma,**kwargs): - return gaussian_filter(I,sigma,**kwargs) - -def imlogfilt(I,sigma,**kwargs): - return -gaussian_laplace(I,sigma,**kwargs) - -def imgradmag(I,sigma): - if len(I.shape) == 2: - dx = imgaussfilt(I,sigma,order=[0,1]) - dy = imgaussfilt(I,sigma,order=[1,0]) - return np.sqrt(dx**2+dy**2) - if len(I.shape) == 3: - dx = imgaussfilt(I,sigma,order=[0,0,1]) - dy = imgaussfilt(I,sigma,order=[0,1,0]) - dz = imgaussfilt(I,sigma,order=[1,0,0]) - return np.sqrt(dx**2+dy**2+dz**2) - -def localstats(I,radius,justfeatnames=False): - ptls = [10,30,50,70,90] - featNames = [] - for i in range(len(ptls)): - featNames.append('locPtl%d' % ptls[i]) - if justfeatnames == True: - return featNames - sI = size(I) - nFeats = len(ptls) - F = np.zeros((sI[0],sI[1],nFeats)) - for i in range(nFeats): - F[:,:,i] = ptlfilt(I,ptls[i],radius) - return F - -def localstats3(I,radius,justfeatnames=False): - ptls = [10,30,50,70,90] - featNames = [] - for i in range(len(ptls)): - featNames.append('locPtl%d' % ptls[i]) - if justfeatnames == True: - return featNames - sI = size(I) - nFeats = len(ptls) - F = np.zeros((sI[0],sI[1],sI[2],nFeats)) - for i in range(nFeats): - F[:,:,:,i] = ptlfilt(I,ptls[i],radius) - return F - -def imderivatives(I,sigmas,justfeatnames=False): - if type(sigmas) is not list: - sigmas = [sigmas] - derivPerSigmaFeatNames = ['d0','dx','dy','dxx','dxy','dyy','normGrad','normHessDiag'] - if justfeatnames == True: - featNames = []; - for i in range(len(sigmas)): - for j in range(len(derivPerSigmaFeatNames)): - featNames.append('derivSigma%d%s' % (sigmas[i],derivPerSigmaFeatNames[j])) - return featNames - nDerivativesPerSigma = len(derivPerSigmaFeatNames) - nDerivatives = len(sigmas)*nDerivativesPerSigma - sI = size(I) - D = np.zeros((sI[0],sI[1],nDerivatives)) - for i in range(len(sigmas)): - sigma = sigmas[i] - dx = imgaussfilt(I,sigma,order=[0,1]) - dy = imgaussfilt(I,sigma,order=[1,0]) - dxx = imgaussfilt(I,sigma,order=[0,2]) - dyy = imgaussfilt(I,sigma,order=[2,0]) - D[:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma) - D[:,:,nDerivativesPerSigma*i+1] = dx - D[:,:,nDerivativesPerSigma*i+2] = dy - D[:,:,nDerivativesPerSigma*i+3] = dxx - D[:,:,nDerivativesPerSigma*i+4] = imgaussfilt(I,sigma,order=[1,1]) - D[:,:,nDerivativesPerSigma*i+5] = dyy - D[:,:,nDerivativesPerSigma*i+6] = np.sqrt(dx**2+dy**2) - D[:,:,nDerivativesPerSigma*i+7] = np.sqrt(dxx**2+dyy**2) - return D - # derivatives are indexed by the last dimension, which is good for ML features but not for visualization, - # in which case the expected dimensions are [plane,channel,y(row),x(col)]; to obtain that ordering, do - # D = np.moveaxis(D,[0,3,1,2],[0,1,2,3]) - -def imderivatives3(I,sigmas,justfeatnames=False): - if type(sigmas) is not list: - sigmas = [sigmas] - - derivPerSigmaFeatNames = ['d0','dx','dy','dz','dxx','dxy','dxz','dyy','dyz','dzz','normGrad','normHessDiag'] - - # derivPerSigmaFeatNames = ['d0','normGrad','normHessDiag'] - - if justfeatnames == True: - featNames = []; - for i in range(len(sigmas)): - for j in range(len(derivPerSigmaFeatNames)): - featNames.append('derivSigma%d%s' % (sigmas[i],derivPerSigmaFeatNames[j])) - return featNames - nDerivativesPerSigma = len(derivPerSigmaFeatNames) - nDerivatives = len(sigmas)*nDerivativesPerSigma - sI = size(I) - D = np.zeros((sI[0],sI[1],sI[2],nDerivatives)) # plane, channel, y, x - for i in range(len(sigmas)): - sigma = sigmas[i] - dx = imgaussfilt(I,sigma,order=[0,0,1]) # z, y, x - dy = imgaussfilt(I,sigma,order=[0,1,0]) - dz = imgaussfilt(I,sigma,order=[1,0,0]) - dxx = imgaussfilt(I,sigma,order=[0,0,2]) - dyy = imgaussfilt(I,sigma,order=[0,2,0]) - dzz = imgaussfilt(I,sigma,order=[2,0,0]) - - D[:,:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma) - D[:,:,:,nDerivativesPerSigma*i+1 ] = dx - D[:,:,:,nDerivativesPerSigma*i+2 ] = dy - D[:,:,:,nDerivativesPerSigma*i+3 ] = dz - D[:,:,:,nDerivativesPerSigma*i+4 ] = dxx - D[:,:,:,nDerivativesPerSigma*i+5 ] = imgaussfilt(I,sigma,order=[0,1,1]) - D[:,:,:,nDerivativesPerSigma*i+6 ] = imgaussfilt(I,sigma,order=[1,0,1]) - D[:,:,:,nDerivativesPerSigma*i+7 ] = dyy - D[:,:,:,nDerivativesPerSigma*i+8 ] = imgaussfilt(I,sigma,order=[1,1,0]) - D[:,:,:,nDerivativesPerSigma*i+9 ] = dzz - D[:,:,:,nDerivativesPerSigma*i+10] = np.sqrt(dx**2+dy**2+dz**2) - D[:,:,:,nDerivativesPerSigma*i+11] = np.sqrt(dxx**2+dyy**2+dzz**2) - - # D[:,:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma) - # D[:,:,:,nDerivativesPerSigma*i+1 ] = np.sqrt(dx**2+dy**2+dz**2) - # D[:,:,:,nDerivativesPerSigma*i+2 ] = np.sqrt(dxx**2+dyy**2+dzz**2) - return D - # derivatives are indexed by the last dimension, which is good for ML features but not for visualization, - # in which case the expected dimensions are [plane,y(row),x(col)]; to obtain that ordering, do - # D = np.moveaxis(D,[2,0,1],[0,1,2]) - -def imfeatures(I=[],sigmaDeriv=1,sigmaLoG=1,locStatsRad=0,justfeatnames=False): - if type(sigmaDeriv) is not list: - sigmaDeriv = [sigmaDeriv] - if type(sigmaLoG) is not list: - sigmaLoG = [sigmaLoG] - derivFeatNames = imderivatives([],sigmaDeriv,justfeatnames=True) - nLoGFeats = len(sigmaLoG) - locStatsFeatNames = [] - if locStatsRad > 1: - locStatsFeatNames = localstats([],locStatsRad,justfeatnames=True) - nLocStatsFeats = len(locStatsFeatNames) - if justfeatnames == True: - featNames = derivFeatNames - for i in range(nLoGFeats): - featNames.append('logSigma%d' % sigmaLoG[i]) - for i in range(nLocStatsFeats): - featNames.append(locStatsFeatNames[i]) - return featNames - nDerivFeats = len(derivFeatNames) - nFeatures = nDerivFeats+nLoGFeats+nLocStatsFeats - sI = size(I) - F = np.zeros((sI[0],sI[1],nFeatures)) - F[:,:,:nDerivFeats] = imderivatives(I,sigmaDeriv) - for i in range(nLoGFeats): - F[:,:,nDerivFeats+i] = imlogfilt(I,sigmaLoG[i]) - if locStatsRad > 1: - F[:,:,nDerivFeats+nLoGFeats:] = localstats(I,locStatsRad) - return F - -def imfeatures3(I=[],sigmaDeriv=2,sigmaLoG=2,locStatsRad=0,justfeatnames=False): - if type(sigmaDeriv) is not list: - sigmaDeriv = [sigmaDeriv] - if type(sigmaLoG) is not list: - sigmaLoG = [sigmaLoG] - derivFeatNames = imderivatives3([],sigmaDeriv,justfeatnames=True) - nLoGFeats = len(sigmaLoG) - locStatsFeatNames = [] - if locStatsRad > 1: - locStatsFeatNames = localstats3([],locStatsRad,justfeatnames=True) - nLocStatsFeats = len(locStatsFeatNames) - if justfeatnames == True: - featNames = derivFeatNames - for i in range(nLoGFeats): - featNames.append('logSigma%d' % sigmaLoG[i]) - for i in range(nLocStatsFeats): - featNames.append(locStatsFeatNames[i]) - return featNames - nDerivFeats = len(derivFeatNames) - nFeatures = nDerivFeats+nLoGFeats+nLocStatsFeats - sI = size(I) - F = np.zeros((sI[0],sI[1],sI[2],nFeatures)) - F[:,:,:,:nDerivFeats] = imderivatives3(I,sigmaDeriv) - for i in range(nLoGFeats): - F[:,:,:,nDerivFeats+i] = imlogfilt(I,sigmaLoG[i]) - if locStatsRad > 1: - F[:,:,:,nDerivFeats+nLoGFeats:] = localstats3(I,locStatsRad) - return F - -def stack2list(S): - L = [] - for i in range(size(S)[2]): - L.append(S[:,:,i]) - return L - -def thrsegment(I,wsBlr,wsThr): # basic threshold segmentation - G = imgaussfilt(I,sigma=(1-wsBlr)+wsBlr*5) # min 1, max 5 - M = G > wsThr - return M \ No newline at end of file