modulecolornorm# Color normalization of microscopy images of Langerhans islets# Jan Kybic, July 2015## written in JuliausingImagesusingImageViewusingColorsusingDebugtypealias ColorType{T}RGB{T}typealias ColorImgType{T,N}Image{ColorType{T},N,Array{ColorType{T},N}}function rgb2yuvimg{T}(img::ColorImgType{T,2})# convert an RGB image into YUV matriximsize=size(img)r=zeros(Float32,(3,imsize[1],imsize[2]))fori=1:length(img)# process all pixelsc=img[i]y=0.299*c.r+0.587*c.g+0.114*c.b# Y - luminosityr[1,i]=yr[2,i]=0.5+0.5*(c.b-y)# u between 0 and 1, not standardr[3,i]=0.5+0.5*(c.r-y)# v between 0 and 1endreturnrendfunction rgb2lbrimg{T}(img::ColorImgType{T,2})# convert an RGB image into luminosity and relative blue and redimsize=size(img)#r=zeros(Float32,(imsize[1],imsize[2]))r=zeros(Float32,(3,imsize[1],imsize[2]))fori=1:length(img)# process all pixelsc=img[i]y=c.r+c.g+c.b# Y - luminosity# J.Schier navrhuje:# r[1,i]=Float32(y)# r[2,i]=Float32(c.b)/y# r[3,i]=Float32(c.r)/yr[1,i]=yr[2,i]=c.b/yr[3,i]=c.r/yendreturnrendfunction yuv2rgb(yuv::Vector)b=2.*yuv[2]-1+yuv[1]r=2.*yuv[3]-1+yuv[1]g=(yuv[1]-0.299*r-0.114*b)/0.587return[r,g,b]endfunction lbr2rgb(l::Vector)b=l[2]*l[1]r=l[3]*l[1]g=l[1]-b-rreturn[r,g,b]endfunction rgb2lbr(x::Vector)l=x[1]+x[2]+x[3]b=x[3]/lr=x[1]/lreturn[l,b,r]endfunction rgb2yuv(rgb::Vector)y=0.299*rgb[1]+0.587*rgb[2]+0.114*rgb[3]# Y - luminosityu=0.5+0.5*(rgb[3]-y)# u between 0 and 1, not standardv=0.5+0.5*(rgb[1]-y)# v between 0 and 1return[y,u,v]endfunction percentile(vec::Vector,perc)# return a threshold value so that "perc" values in "vec" are smaller and 1-perc are larger@assert(perc>=0.&&perc<=1.)sv=sort(vec)ind=min(round(Int,floor(length(vec)*perc))+1,length(vec))returnsv[ind]endfunction hist2(m::Matrix,nbins)# create a 2D histogram from a (2,n) matrix, assume the values between 0 and 1 with nbins in each dimensionn=size(m,2)@assert(size(m,1)==2)h=zeros(UInt32,nbins,nbins)nbinsf=float(nbins)fori=1:niu=min(round(Int,floor(m[1,i]*nbinsf))+1,nbins)iv=min(round(Int,floor(m[2,i]*nbinsf))+1,nbins)h[iu,iv]+=1endreturnhendfunction hist2max(h::Matrix)# find the maximum element of the histogram and return the corresponding valuesnbins=size(h,1)rnbins=1./float(nbins)@assert(nbins==size(h,2))(iu,iv)=ind2sub(size(h),indmax(h))u=(iu-0.5)*rnbinsv=(iv-0.5)*rnbinsreturn(u,v)endfunction rgbimgmax{T}(img::ColorImgType{T,2})# return the maximum value in all channels of a RGB imagereturnmaximum(reshape(separate(img),size(img,1)*size(img,2),3),1).data'endfunction rgbimgpercentile{T}(img::ColorImgType{T,2};perc=0.01)# almost like rgbimg but ignore the highest 'perc' pixels (relative number)q=reshape(separate(img),size(img,1)*size(img,2),3).datareturn[percentile(q[:,1],1-perc),percentile(q[:,2],1-perc),percentile(q[:,3],1-perc)]endfunction piecewiselintransf(x0,y0,x)# calculates y=f(x), so that f(0)=0, f(1)=1, f(x0)=y0 and between these values the function is linearifx>x0return(x-x0)/(1-x0)*(1-y0)+y0elsereturnx/x0*y0endend@debugfunction inventbluechannel{T}(img::ColorImgType{T,2};minval=0.2)# if the maximum blue channel value is smaller than minval, invent the blue value#maxrgb=rgbimgmax(img)maxrgb=rgbimgpercentile(img,perc=0.1)#println("maxrgb=$maxrgb")ifmaxrgb[3]<minvalprintln("Inventing blue channel values")r=similar(img)# output imagefori=1:length(img)r[i]=RGB(img[i].r,img[i].g,img[i].g)endreturnrendreturnimgend@debugfunction normalizeimg{T,N}(img::ColorImgType{T,N};topvals=0.1,nbins=100,method="max",corrmethod="ignoreoutliers",inventblue=true)# take an RGB image and change its colors so that the background is white# assume that most pixels in the image have the background pixels and that the "topvals" brightest pixels# belong to the background. "nbins" is the number of bins of a histogram# topvals is between 0 and 1.# method can be 'max' (most frequent color),# 'mean' (mean color in relative RGB space)# 'meanRGB' (mean color in RGB space)# corrmethod can be 'norm',# or 'hueonly'# or 'limit'# or 'clip'@assert(length(size(img))==2)ifinventblueimg=inventbluechannel(img)endr=similar(img)# output imagelbr=rgb2lbrimg(img)# convert to luminosity and relative blue and red components# only consider topvals*length(img) brightest valueslvec=reshape(lbr[1,:,:],length(img))threshold=percentile(lvec,topvals)inds=lvec.>threshold# indices to be consideredbr=reshape(lbr[2:3,:,:],2,length(img))br=br[:,inds]# BR values of the brightest pixelsifmethod=="max"h=hist2(br,nbins)# calculate histogram and find its maximumbb,br=hist2max(h)# coordinates of the maximum = background color BR componentsbrgb=lbr2rgb([1.,bb,br])# background color in RGBelseifmethod=="mean"# use the meanbbr=mean(br,2)bb=bbr[1];br=bbr[2]brgb=lbr2rgb([1.,bb,br])# background color in RGBelseifmethod=="meanRGB"# use the mean calculated in RGBrgbv=reshape(float32(separate(img).data),size(img,1)*size(img,2),3)brgb=mean(rgbv[inds,:],1)elseerror("Unknown method `$(method)'")end# now the bacground color is in brgb# calculate gain in RGB channels so that the background color "brgb" becomes whiteg=1./float(brgb)ifcorrmethod=="norm"# adjust the gain so that the maximum value is onemaxrgb=rgbimgmax(img)g*=1./maximum(maxrgb.*g)fori=1:length(img)# process all pixelsr[i]=RGB(g[1]*img[i].r,g[2]*img[i].g,g[3]*img[i].b)endelseifcorrmethod=="ignoreoutliers"# adjust the gain so that the maximum value is one except at very few placesmaxrgb=rgbimgpercentile(img,perc=0.01)g*=1./maximum(maxrgb.*g)fori=1:length(img)# process all pixelsr[i]=RGB(min(g[1]*img[i].r,1.0),min(g[2]*img[i].g,1.0),min(g[3]*img[i].b,1.0))endelseifcorrmethod=="hueonly"# adjust the gain so that the sum is onefori=1:length(img)# process all pixelsy=img[i].r+img[i].g+img[i].brn=g[1]*img[i].r;gn=g[2]*img[i].g;bn=g[3]*img[i].b;# scaling factor - keep intensity if possible but do not clipc=min(y/(rn+gn+bn),1./rn,1./gn,1./bn)r[i]=RGB(c*rn,c*gn,c*bn)endelseifcorrmethod=="hueonlyclip"# adjust the gain so that the sum is onefori=1:length(img)# process all pixelsy=img[i].r+img[i].g+img[i].brn=g[1]*img[i].r;gn=g[2]*img[i].g;bn=g[3]*img[i].b;# scaling factor - keep intensity if possible and clip if not possiblec=y/(rn+gn+bn)r[i]=RGB(min(1.,c*rn),min(1.,c*gn),min(1.,c*bn))endelseifcorrmethod=="clip"fori=1:length(img)# process all pixelsr[i]=RGB(min(g[1]*img[i].r,1.0),min(g[2]*img[i].g,1.0),min(g[3]*img[i].b,1.0))endelseifcorrmethod=="piecewise"t=sum(brgb)/3.# target grayfori=1:length(img)# process all pixelsr[i]=RGB(piecewiselintransf(brgb[1],t,img[i].r),piecewiselintransf(brgb[2],t,img[i].g),piecewiselintransf(brgb[3],t,img[i].b))endelseerror("Unknown corrmethod `$(corrmethod)'")endreturnrendfunction test_normalizeimg()# try color normalization on one imagea=imread(expanduser("~/work/IKEM/good_clinical/iz130305btxN1-10x.png"))b=normalizeimg(a,)ImageView.view(b)endfunction normalize_from_directory(inpdir,outdir)# normalize all images in a directory inpdir and save them to outpdirassert(isdir(inpdir))assert(isdir(outdir))files=readdir(inpdir)picture_extensions=Set([".png",".jpg",".bmp"])forfilenameinfilesfn,ext=splitext(filename)# process only image filesifin(ext,picture_extensions)infilename=joinpath(inpdir,filename)outfilename=joinpath(outdir,filename)println("Processing file `$(infilename)' to `$(outfilename)'.")ifisfile(outfilename)||isdir(outfilename)println("File exists, skipping.")continueenda=imread(infilename)println("Read")b=atryb=normalizeimg(a)println("Processed")catchprintln("Processing failed, reverting.")endimwrite(b,outfilename)println("Written")elseprintln("Ignoring file `$(file)'.")endendendfunction test_colornorm()# normalize all images from a given directory# can be launched e.g. by import colornorm ; colornorm.test_colornorm()normalize_from_directory(expanduser("~/work/IKEM/good_clinical/"),expanduser("~/work/IKEM/good_clinical_colornorm"))endfunction normalize_training_set()normalize_from_directory("/mnt/tmp2/microscopy/IKEM/GT_IMAGE_TRENOVACI_SADA/images",expanduser("~/work/IKEM/images_colornorm"))endfunction normalize_all_clinical()normalize_from_directory("/mnt/tmp2/microscopy/IKEM/data20150730_all_clinical",expanduser("~/work/IKEM/data20150730_all_clinical_colornorm"))endfunction normalize_independent_donors()srcprefix="/mnt/tmp/microscopy/IKEM/Independent Donors2015"dstprefix="/mnt/tmp/microscopy/IKEM/Independent Donors2015_colornorm"forsuffixin["Donor7_iz150821d0","Donor8_iz150724d2","Donor9_iz150423d0"]normalize_from_directory(joinpath(srcprefix,suffix),joinpath(dstprefix,suffix))endendfunction normalize_independent_donors2()srcprefix="/mnt/tmp/microscopy/IKEM/Independent_Donors2015_2"dstprefix="/mnt/tmp/microscopy/IKEM/Independent_Donors2015_2_colornorm"normalize_from_directory(srcprefix,dstprefix)endfunction normalize_GT_IMAGE_PRO_BAREV_NORM()normalize_from_directory(expanduser("~/work/IKEM/GT_IMAGE_PRO_BAREV_NORM"),expanduser("~/work/IKEM/GT_IMAGE_PRO_BAREV_NORM_colornorm"))endfunction normalize_151015()normalize_from_directory("/mnt/tmp2/microscopy/IKEM/nova_izolace_15_10_15","/mnt/tmp2/microscopy/IKEM/nova_izolace_15_10_15_colornorm")endfunction normalize_241015()normalize_from_directory("/mnt/tmp2/microscopy/IKEM/nova_izolace_24_10_15","/mnt/tmp2/microscopy/IKEM/nova_izolace_24_10_15_colornorm")endfunction normalize_011115()normalize_from_directory("/mnt/tmp2/microscopy/IKEM/nova_izolace_1_11_15","/mnt/tmp2/microscopy/IKEM/nova_izolace_1_11_15_colornorm")normalize_from_directory("/mnt/tmp2/microscopy/IKEM/PURITY_1_11_2015","/mnt/tmp2/microscopy/IKEM/PURITY_1_11_2015_colornorm")endfunction normalize_thisdir()normalize_from_directory("inpimgs","outpimgs")endend