{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\n\nprint(os.listdir(\"../input\"))\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import os, sys\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport skimage.io\nfrom skimage.transform import resize\nfrom imgaug import augmenters as iaa\nfrom tqdm import tqdm\nimport PIL\nfrom PIL import Image, ImageOps\nimport cv2\nfrom sklearn.utils import class_weight, shuffle\nfrom keras.losses import binary_crossentropy\nfrom keras.applications.resnet50 import preprocess_input\nimport keras.backend as K\nimport tensorflow as tf\nfrom sklearn.metrics import f1_score, fbeta_score\nfrom keras.utils import Sequence\nfrom keras.utils import to_categorical\nfrom sklearn.model_selection import train_test_split\n\nWORKERS = 2\nCHANNEL = 3\n\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nIMG_SIZE = 512\nNUM_CLASSES = 5\nSEED = 77\nTRAIN_NUM = 1000 # use 1000 when you just want to explore new idea, use -1 for full train","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_df = pd.read_csv('../input/aptos2019-blindness-detection/train.csv')\ntest_df = pd.read_csv('../input/aptos2019-blindness-detection/test.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import cmath as cm\nimport numpy as np","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class circularHOGExtractor():\n    \"\"\"\n    This method takes in a single image and extracts rotation invariant HOG features\n    following the approach in this paper: \n    Liu, Kun, et al. \"Rotation-invariant HOG descriptors using fourier analysis in polar and spherical coordinates.\"\n    International Journal of Computer Vision 106.3 (2014): 342-364.\n    \"\"\"\n    def __init__(self, bins=4, size=6, max_freq=4):\n\n        # number of bins in the radial direction for large scale features\n        self.mNBins = bins\n        # size of bin in pixels, this sets the required radius for the image = bins*size\n        self.mNSize = size\n        # number of fourier modes that will be used (0:modes-1)\n        self.mNMaxFreq = max_freq \n\n        mf = self.mNMaxFreq+1\n        self.mNCount = 2*(bins-1) * (mf + 2*(np.dot([mf - i for i in range(mf)] , range(mf))  ))\n        # create a list to store kernels for regional descriptors based on circular harmonics\n        self.ciKernel = []\n\n        # first create the central region \n        [x,y]=np.meshgrid(range(-self.mNSize+1,self.mNSize),range(-self.mNSize+1,self.mNSize))\n        z = x + 1j*y\n        kernel = self.mNSize - np.abs(z)\n        kernel[kernel < 0] = 0\n        kernel = kernel/sum(sum(kernel))\n\n #       self.ciKernel.append(kernel)\n\n        # next build the internal regions - (bins-1) concentric circles\n        modes = range(0, self.mNMaxFreq+1)\n        scale = range(2, self.mNBins+1)\n\n        for s in scale:\n            r = int(self.mNSize * s)\n            ll = range(1-r,r)\n            [x,y] = np.meshgrid(ll,ll)\n            z = x + 1j*y\n            phase_z = np.angle(z);\n                            \n            for k in modes:\n                kernel = self.mNSize - np.abs(np.abs(z) - (r-self.mNSize)) \n                kernel[kernel < 0] = 0\n                kernel = np.multiply(kernel,np.exp(1j*phase_z*k))\n                sa = np.ravel(np.abs(kernel))\n                kernel = kernel / np.sqrt(np.sum(np.multiply(sa,sa)))\n\n                self.ciKernel.append(kernel)\n\n\n\n    def extract(self, img):\n        I = img.astype(float)/255.0\n#      I = (I-I.mean())/I.std()\n\n        # size and centre of image\n        (nx, ny) = I.shape\n        cx = int(round(0.5*nx))\n        cy = int(round(0.5*ny))\n\n        # compute gradient with a central difference method and store in complex form\n        (dy, dx) = np.gradient(I)\n        dz = dx + 1j*dy\n\n        # compute magnitude/phase of complex numbers\n        phi = np.angle(dz)\n        r = np.abs(dz)\n#       r = r/(r.std()+0.0001)\n\n\n        # create an empty array for storing the dfft of the orientation vector\n        histF = np.zeros([nx, ny, self.mNMaxFreq+1])+0j\n\n        # take the dfft of the orientation vector up to order MaxFreq\n        # positive values of k only since negative values give conjugate\n        for k in range(0,self.mNMaxFreq+1):\n            histF[:,:,k] = np.multiply(np.exp( -1j * (k) * phi) , r+0j)\n        \n\n        # compute regional descriptors by convolutions (these descriptors are not rotation invariant)\n        fHOG = np.zeros([self.mNCount])\n        scale = range(0, self.mNBins-1)\n        f_index = 0\n        for s in scale:\n            allVals = np.zeros((self.mNMaxFreq+1,self.mNMaxFreq+1),dtype=np.complex64)\n            for freq in range(0,self.mNMaxFreq+1):\n                template = self.ciKernel[s*(self.mNMaxFreq+1)+freq]\n                (tnx, tny) = template.shape\n                tnx2 = int(round(0.5*tnx))\n                for k in range(0,self.mNMaxFreq+1):\n                    allVals[freq,k] = np.sum(np.sum(np.multiply(histF[cx-tnx2:cx-tnx2+tnx,cy-tnx2:cy-tnx2+tnx,k],template)))\n            for (x,y), val in np.ndenumerate(allVals):\n                if x==y:\n                    fHOG[f_index]=val.real\n                    f_index+=1\n                    fHOG[f_index]=val.imag\n                    f_index+=1\n                else:\n                    for (x1,y1), val1 in np.ndenumerate(allVals):\n                        if x1<x: continue\n                        if y1<y: continue\n                        if (x-y)==(x1-y1):\n                            fHOG[f_index]=(val*val1.conjugate()).real\n                            f_index+=1\n                            fHOG[f_index]=(val*val1.conjugate()).imag\n                            f_index+=1\n\n        return fHOG.tolist()\n\n\n    def prepareExtract(self, img):\n        I = img.astype(float)/255.0\n#      I = (I-I.mean())/I.std()\n\n        # size and centre of image\n        (nx, ny) = I.shape\n        # compute gradient with a central difference method and store in complex form\n        (dy, dx) = np.gradient(I)\n        dz = dx + 1j*dy\n\n        # compute magnitude/phase of complex numbers\n        phi = np.angle(dz)\n        r = np.abs(dz)\n #       r = r/(r.mean()+0.001)\n\n\n        # create an empty array for storing the dfft of the orientation vector\n        histF = np.zeros([nx, ny, self.mNMaxFreq+1])+0j\n\n        # take the dfft of the orientation vector up to order MaxFreq\n        # positive values of k only since negative values give conjugate\n        for k in range(0,self.mNMaxFreq+1):\n            histF[:,:,k] = np.multiply(np.exp( -1j * (k) * phi) , r+0j)\n\n        return histF\n        \n    def denseExtract(self, histF, positions, N):\n #       I = img.astype(float)/255.0\n#      I = (I-I.mean())/I.std()\n\n        # size and centre of image\n        (nx, ny, kk) = histF.shape\n        \n        features = np.zeros((N,self.mNCount),dtype=np.float32)\n        scale = range(0, self.mNBins-1)\n        for p in range(N):\n            cx = positions[p,0]+1\n            cy = positions[p,1]+1\n            if cx<self.mNBins*self.mNSize: continue\n            if cy<self.mNBins*self.mNSize: continue\n            if cx> nx - self.mNBins*self.mNSize: continue\n            if cy> ny - self.mNBins*self.mNSize: continue\n            \n            f_index = 0\n            for s in scale:\n                allVals = np.zeros((self.mNMaxFreq+1,self.mNMaxFreq+1),dtype=np.complex64)\n                \n                for freq in range(0,self.mNMaxFreq+1):\n                    template = self.ciKernel[s*(self.mNMaxFreq+1)+freq]\n                    (tnx, tny) = template.shape\n                    tnx2 = int(round(0.5*tnx))\n                    for k in range(0,self.mNMaxFreq+1):\n                        allVals[freq,k] = np.sum(np.sum(np.multiply(histF[cx-tnx2:cx-tnx2+tnx,cy-tnx2:cy-tnx2+tnx,k],template)))\n                        #if p==2193 and freq==0 and s==0:\n                        #        print k\n                        #        for kk in histF[cx-tnx2:cx-tnx2+tnx,cy-tnx2:cy-tnx2+tnx,k]:\n                        #            for jj in kk:\n                        #                print jj.real\n                \n                \n                for (x,y), val in np.ndenumerate(allVals):\n                    if x==y:\n                        features[p,f_index]=val.real\n                        f_index+=1\n                        features[p,f_index]=val.imag\n                        f_index+=1\n\n                    else:\n                        for (x1,y1), val1 in np.ndenumerate(allVals):\n                            if x1<x: continue\n                            if y1<y: continue\n                            if (x-y)==(x1-y1):\n                                features[p,f_index]=(val*val1.conjugate()).real\n                                f_index+=1\n                                features[p,f_index]=(val*val1.conjugate()).imag\n                                f_index+=1\n\n        \n        return features\n\n#        print \"diff to original array:\"\n#        print features[0], fHOG[0]\n#        print np.max(np.abs(features-fHOG))\n\n        return fHOG.tolist()\n\n    \n    def getFieldNames(self):\n        \"\"\"\n        Return the names of all of the length and angle fields. \n        \"\"\"\n        retVal = []\n        for i in range(0,self.mNCount):\n            name = \"Length\"+str(i)\n            retVal.append(name)\n                        \n        return retVal\n        \"\"\"\n        This method gives the names of each field in the feature vector in the\n        order in which they are returned. For example, 'xpos' or 'width'\n        \"\"\"\n\n    def getNumFields(self):\n        \"\"\"\n        This method returns the total number of fields in the feature vector.\n        \"\"\"\n        return self.mNCount\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df0 = train_df[train_df['diagnosis']==0]\ndf1 = train_df[train_df['diagnosis']==1]\ndf2 = train_df[train_df['diagnosis']==2]\ndf3 = train_df[train_df['diagnosis']==3]\ndf4 = train_df[train_df['diagnosis']==4]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df0.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dflist = [df0,df1,df2,df3,df4]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dflist[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class circularHOGExtractor():\n    \"\"\"\n    This method takes in a single image and extracts rotation invariant HOG features\n    following the approach in this paper: \n    Liu, Kun, et al. \"Rotation-invariant HOG descriptors using fourier analysis in polar and spherical coordinates.\"\n    International Journal of Computer Vision 106.3 (2014): 342-364.\n    \"\"\"\n    def __init__(self, bins=4, size=6, max_freq=4):\n\n        # number of bins in the radial direction for large scale features\n        self.mNBins = bins\n        # size of bin in pixels, this sets the required radius for the image = bins*size\n        self.mNSize = size\n        # number of fourier modes that will be used (0:modes-1)\n        self.mNMaxFreq = max_freq \n\n        mf = self.mNMaxFreq+1\n        self.mNCount = 2*(bins-1) * (mf + 2*(np.dot([mf - i for i in range(mf)] , range(mf))  ))\n        # create a list to store kernels for regional descriptors based on circular harmonics\n        self.ciKernel = []\n\n        # first create the central region \n        [x,y]=np.meshgrid(range(-self.mNSize+1,self.mNSize),range(-self.mNSize+1,self.mNSize))\n        z = x + 1j*y\n        kernel = self.mNSize - np.abs(z)\n        kernel[kernel < 0] = 0\n        kernel = kernel/sum(sum(kernel))\n\n #       self.ciKernel.append(kernel)\n\n        # next build the internal regions - (bins-1) concentric circles\n        modes = range(0, self.mNMaxFreq+1)\n        scale = range(2, self.mNBins+1)\n\n        for s in scale:\n            r = int(self.mNSize * s)\n            ll = range(1-r,r)\n            [x,y] = np.meshgrid(ll,ll)\n            z = x + 1j*y\n            phase_z = np.angle(z);\n                            \n            for k in modes:\n                kernel = self.mNSize - np.abs(np.abs(z) - (r-self.mNSize)) \n                kernel[kernel < 0] = 0\n                kernel = np.multiply(kernel,np.exp(1j*phase_z*k))\n                sa = np.ravel(np.abs(kernel))\n                kernel = kernel / np.sqrt(np.sum(np.multiply(sa,sa)))\n\n                self.ciKernel.append(kernel)\n\n\n\n    def extract(self, img):\n        I = img.astype(float)/255.0\n#      I = (I-I.mean())/I.std()\n\n        # size and centre of image\n        (nx, ny) = I.shape\n        cx = int(round(0.5*nx))\n        cy = int(round(0.5*ny))\n\n        # compute gradient with a central difference method and store in complex form\n        (dy, dx) = np.gradient(I)\n        dz = dx + 1j*dy\n\n        # compute magnitude/phase of complex numbers\n        phi = np.angle(dz)\n        r = np.abs(dz)\n#       r = r/(r.std()+0.0001)\n\n\n        # create an empty array for storing the dfft of the orientation vector\n        histF = np.zeros([nx, ny, self.mNMaxFreq+1])+0j\n\n        # take the dfft of the orientation vector up to order MaxFreq\n        # positive values of k only since negative values give conjugate\n        for k in range(0,self.mNMaxFreq+1):\n            histF[:,:,k] = np.multiply(np.exp( -1j * (k) * phi) , r+0j)\n        \n\n        # compute regional descriptors by convolutions (these descriptors are not rotation invariant)\n        fHOG = np.zeros([self.mNCount])\n        scale = range(0, self.mNBins-1)\n        f_index = 0\n        for s in scale:\n            allVals = np.zeros((self.mNMaxFreq+1,self.mNMaxFreq+1),dtype=np.complex64)\n            for freq in range(0,self.mNMaxFreq+1):\n                template = self.ciKernel[s*(self.mNMaxFreq+1)+freq]\n                (tnx, tny) = template.shape\n                tnx2 = int(round(0.5*tnx))\n                for k in range(0,self.mNMaxFreq+1):\n                    allVals[freq,k] = np.sum(np.sum(np.multiply(histF[cx-tnx2:cx-tnx2+tnx,cy-tnx2:cy-tnx2+tnx,k],template)))\n            for (x,y), val in np.ndenumerate(allVals):\n                if x==y:\n                    fHOG[f_index]=val.real\n                    f_index+=1\n                    fHOG[f_index]=val.imag\n                    f_index+=1\n                else:\n                    for (x1,y1), val1 in np.ndenumerate(allVals):\n                        if x1<x: continue\n                        if y1<y: continue\n                        if (x-y)==(x1-y1):\n                            fHOG[f_index]=(val*val1.conjugate()).real\n                            f_index+=1\n                            fHOG[f_index]=(val*val1.conjugate()).imag\n                            f_index+=1\n\n        return fHOG.tolist()\n\n\n    def prepareExtract(self, img):\n        I = img.astype(float)/255.0\n#      I = (I-I.mean())/I.std()\n\n        # size and centre of image\n        (nx, ny) = I.shape\n        # compute gradient with a central difference method and store in complex form\n        (dy, dx) = np.gradient(I)\n        dz = dx + 1j*dy\n\n        # compute magnitude/phase of complex numbers\n        phi = np.angle(dz)\n        r = np.abs(dz)\n #       r = r/(r.mean()+0.001)\n\n\n        # create an empty array for storing the dfft of the orientation vector\n        histF = np.zeros([nx, ny, self.mNMaxFreq+1])+0j\n\n        # take the dfft of the orientation vector up to order MaxFreq\n        # positive values of k only since negative values give conjugate\n        for k in range(0,self.mNMaxFreq+1):\n            histF[:,:,k] = np.multiply(np.exp( -1j * (k) * phi) , r+0j)\n\n        return histF\n        \n    def denseExtract(self, histF, positions, N):\n #       I = img.astype(float)/255.0\n#      I = (I-I.mean())/I.std()\n\n        # size and centre of image\n        (nx, ny, kk) = histF.shape\n        \n        features = np.zeros((N,self.mNCount),dtype=np.float32)\n        scale = range(0, self.mNBins-1)\n        for p in range(N):\n            cx = positions[p,0]+1\n            cy = positions[p,1]+1\n            if cx<self.mNBins*self.mNSize: continue\n            if cy<self.mNBins*self.mNSize: continue\n            if cx> nx - self.mNBins*self.mNSize: continue\n            if cy> ny - self.mNBins*self.mNSize: continue\n            \n            f_index = 0\n            for s in scale:\n                allVals = np.zeros((self.mNMaxFreq+1,self.mNMaxFreq+1),dtype=np.complex64)\n                \n                for freq in range(0,self.mNMaxFreq+1):\n                    template = self.ciKernel[s*(self.mNMaxFreq+1)+freq]\n                    (tnx, tny) = template.shape\n                    tnx2 = int(round(0.5*tnx))\n                    for k in range(0,self.mNMaxFreq+1):\n                        allVals[freq,k] = np.sum(np.sum(np.multiply(histF[cx-tnx2:cx-tnx2+tnx,cy-tnx2:cy-tnx2+tnx,k],template)))\n                        #if p==2193 and freq==0 and s==0:\n                        #        print k\n                        #        for kk in histF[cx-tnx2:cx-tnx2+tnx,cy-tnx2:cy-tnx2+tnx,k]:\n                        #            for jj in kk:\n                        #                print jj.real\n                \n                \n                for (x,y), val in np.ndenumerate(allVals):\n                    if x==y:\n                        features[p,f_index]=val.real\n                        f_index+=1\n                        features[p,f_index]=val.imag\n                        f_index+=1\n\n                    else:\n                        for (x1,y1), val1 in np.ndenumerate(allVals):\n                            if x1<x: continue\n                            if y1<y: continue\n                            if (x-y)==(x1-y1):\n                                features[p,f_index]=(val*val1.conjugate()).real\n                                f_index+=1\n                                features[p,f_index]=(val*val1.conjugate()).imag\n                                f_index+=1\n\n        \n        return features\n\n#        print \"diff to original array:\"\n#        print features[0], fHOG[0]\n#        print np.max(np.abs(features-fHOG))\n\n        return fHOG.tolist()\n\n    \n    def getFieldNames(self):\n        \"\"\"\n        Return the names of all of the length and angle fields. \n        \"\"\"\n        retVal = []\n        for i in range(0,self.mNCount):\n            name = \"Length\"+str(i)\n            retVal.append(name)\n                        \n        return retVal\n        \"\"\"\n        This method gives the names of each field in the feature vector in the\n        order in which they are returned. For example, 'xpos' or 'width'\n        \"\"\"\n\n    def getNumFields(self):\n        \"\"\"\n        This method returns the total number of fields in the feature vector.\n        \"\"\"\n        return self.mNCount\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.feature_extraction.image import extract_patches_2d","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(df3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"samplesize = 100\ndef extract_hog(df):\n    featurelist = []\n    \n    for counter,j in enumerate(df['id_code'].values):\n        if counter%10 == 0:\n            print(counter)\n        img = cv2.imread('../input/aptos2019-blindness-detection/train_images/'+ j + '.png')\n        grayimg = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)\n        grayimg = cv2.resize(grayimg, (IMG_SIZE,IMG_SIZE))\n        patches = extract_patches_2d(grayimg, (64,64))\n        numbers = [i for i in range(len(patches))]\n        randomindex = np.random.choice(numbers, samplesize)\n        sample_patches = patches[randomindex]\n        for i in sample_patches:\n            extractor = circularHOGExtractor()\n            hogfeature = extractor.extract(i)\n            featurelist.append(hogfeature)    \n        #cv2.imread('../input/aptos2019-blindness-detection/train_images/'+ j + '.png').shape\n    return featurelist","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"featurelist = extract_hog(df3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"featurearr = np.array(featurelist)\nfeaturearr.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.cluster import KMeans","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cluster_num =50\nkm = KMeans(n_clusters=cluster_num)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"clusterlist = []\nfor i in np.arange(0,len(featurelist),samplesize):\n    km.fit(featurearr[i:i+samplesize])\n    tmp = km.predict(featurearr[i:i+samplesize])\n    clusterlist.append(tmp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"clusterlist[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"hist = []\nfor i in clusterlist:\n    tmp = [np.sum(i==j) for j in range(cluster_num)]\n    hist.append(tmp)    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(hist)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"len(hist[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#plt.figure(figsize=(30,50))\nfor i in np.random.choice([j for j in range(len(hist))],10):\n    plt.bar([i for i in range(cluster_num)],hist[i])\n    #plt.hist([i for i in range(cluster_num)],hist[i])\n    #sns.distplot(hist[i])\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}