diff toolbox/PartitionOfImage.py @ 0:99308601eaa6 draft

"planemo upload for repository https://github.com/ohsu-comp-bio/UNetCoreograph commit fb90660a1805b3f68fcff80d525b5459c3f7dfd6-dirty"

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/PartitionOfImage.py	Wed May 19 21:34:38 2021 +0000
@@ -0,0 +1,305 @@
+import numpy as np
+from toolbox.imtools import *
+# from toolbox.ftools import *
+# import sys
+
+class PI2D:
+    Image = None
+    PaddedImage = None
+    PatchSize = 128
+    Margin = 14
+    SubPatchSize = 100
+    PC = None # patch coordinates
+    NumPatches = 0
+    Output = None
+    Count = None
+    NR = None
+    NC = None
+    NRPI = None
+    NCPI = None
+    Mode = None
+    W = None
+
+    def setup(image,patchSize,margin,mode):
+        PI2D.Image = image
+        PI2D.PatchSize = patchSize
+        PI2D.Margin = margin
+        subPatchSize = patchSize-2*margin
+        PI2D.SubPatchSize = subPatchSize
+
+        W = np.ones((patchSize,patchSize))
+        W[[0,-1],:] = 0
+        W[:,[0,-1]] = 0
+        for i in range(1,2*margin):
+            v = i/(2*margin)
+            W[i,i:-i] = v
+            W[-i-1,i:-i] = v
+            W[i:-i,i] = v
+            W[i:-i,-i-1] = v
+        PI2D.W = W
+
+        if len(image.shape) == 2:
+            nr,nc = image.shape
+        elif len(image.shape) == 3: # multi-channel image
+            nz,nr,nc = image.shape
+
+        PI2D.NR = nr
+        PI2D.NC = nc
+
+        npr = int(np.ceil(nr/subPatchSize)) # number of patch rows
+        npc = int(np.ceil(nc/subPatchSize)) # number of patch cols
+
+        nrpi = npr*subPatchSize+2*margin # number of rows in padded image 
+        ncpi = npc*subPatchSize+2*margin # number of cols in padded image 
+
+        PI2D.NRPI = nrpi
+        PI2D.NCPI = ncpi
+
+        if len(image.shape) == 2:
+            PI2D.PaddedImage = np.zeros((nrpi,ncpi))
+            PI2D.PaddedImage[margin:margin+nr,margin:margin+nc] = image
+        elif len(image.shape) == 3:
+            PI2D.PaddedImage = np.zeros((nz,nrpi,ncpi))
+            PI2D.PaddedImage[:,margin:margin+nr,margin:margin+nc] = image
+
+        PI2D.PC = [] # patch coordinates [r0,r1,c0,c1]
+        for i in range(npr):
+            r0 = i*subPatchSize
+            r1 = r0+patchSize
+            for j in range(npc):
+                c0 = j*subPatchSize
+                c1 = c0+patchSize
+                PI2D.PC.append([r0,r1,c0,c1])
+
+        PI2D.NumPatches = len(PI2D.PC)
+        PI2D.Mode = mode # 'replace' or 'accumulate'
+
+    def getPatch(i):
+        r0,r1,c0,c1 = PI2D.PC[i]
+        if len(PI2D.PaddedImage.shape) == 2:
+            return PI2D.PaddedImage[r0:r1,c0:c1]
+        if len(PI2D.PaddedImage.shape) == 3:
+            return PI2D.PaddedImage[:,r0:r1,c0:c1]
+
+    def createOutput(nChannels):
+        if nChannels == 1:
+            PI2D.Output = np.zeros((PI2D.NRPI,PI2D.NCPI),np.float16)
+        else:
+            PI2D.Output = np.zeros((nChannels,PI2D.NRPI,PI2D.NCPI),np.float16)
+        if PI2D.Mode == 'accumulate':
+            PI2D.Count = np.zeros((PI2D.NRPI,PI2D.NCPI),np.float16)
+
+    def patchOutput(i,P):
+        r0,r1,c0,c1 = PI2D.PC[i]
+        if PI2D.Mode == 'accumulate':
+            PI2D.Count[r0:r1,c0:c1] += PI2D.W
+        if len(P.shape) == 2:
+            if PI2D.Mode == 'accumulate':
+                PI2D.Output[r0:r1,c0:c1] += np.multiply(P,PI2D.W)
+            elif PI2D.Mode == 'replace':
+                PI2D.Output[r0:r1,c0:c1] = P
+        elif len(P.shape) == 3:
+            if PI2D.Mode == 'accumulate':
+                for i in range(P.shape[0]):
+                    PI2D.Output[i,r0:r1,c0:c1] += np.multiply(P[i,:,:],PI2D.W)
+            elif PI2D.Mode == 'replace':
+                PI2D.Output[:,r0:r1,c0:c1] = P
+
+    def getValidOutput():
+        margin = PI2D.Margin
+        nr, nc = PI2D.NR, PI2D.NC
+        if PI2D.Mode == 'accumulate':
+            C = PI2D.Count[margin:margin+nr,margin:margin+nc]
+        if len(PI2D.Output.shape) == 2:
+            if PI2D.Mode == 'accumulate':
+                return np.divide(PI2D.Output[margin:margin+nr,margin:margin+nc],C)
+            if PI2D.Mode == 'replace':
+                return PI2D.Output[margin:margin+nr,margin:margin+nc]
+        if len(PI2D.Output.shape) == 3:
+            if PI2D.Mode == 'accumulate':
+                for i in range(PI2D.Output.shape[0]):
+                    PI2D.Output[i,margin:margin+nr,margin:margin+nc] = np.divide(PI2D.Output[i,margin:margin+nr,margin:margin+nc],C)
+            return PI2D.Output[:,margin:margin+nr,margin:margin+nc]
+
+
+    def demo():
+        I = np.random.rand(128,128)
+        # PI2D.setup(I,128,14)
+        PI2D.setup(I,64,4,'replace')
+
+        nChannels = 2
+        PI2D.createOutput(nChannels)
+
+        for i in range(PI2D.NumPatches):
+            P = PI2D.getPatch(i)
+            Q = np.zeros((nChannels,P.shape[0],P.shape[1]))
+            for j in range(nChannels):
+                Q[j,:,:] = P
+            PI2D.patchOutput(i,Q)
+
+        J = PI2D.getValidOutput()
+        J = J[0,:,:]
+
+        D = np.abs(I-J)
+        print(np.max(D))
+
+        K = cat(1,cat(1,I,J),D)
+        imshow(K)
+
+
+class PI3D:
+    Image = None
+    PaddedImage = None
+    PatchSize = 128
+    Margin = 14
+    SubPatchSize = 100
+    PC = None # patch coordinates
+    NumPatches = 0
+    Output = None
+    Count = None
+    NR = None # rows
+    NC = None # cols
+    NZ = None # planes
+    NRPI = None
+    NCPI = None
+    NZPI = None
+    Mode = None
+    W = None
+
+    def setup(image,patchSize,margin,mode):
+        PI3D.Image = image
+        PI3D.PatchSize = patchSize
+        PI3D.Margin = margin
+        subPatchSize = patchSize-2*margin
+        PI3D.SubPatchSize = subPatchSize
+
+        W = np.ones((patchSize,patchSize,patchSize))
+        W[[0,-1],:,:] = 0
+        W[:,[0,-1],:] = 0
+        W[:,:,[0,-1]] = 0
+        for i in range(1,2*margin):
+            v = i/(2*margin)
+            W[[i,-i-1],i:-i,i:-i] = v
+            W[i:-i,[i,-i-1],i:-i] = v
+            W[i:-i,i:-i,[i,-i-1]] = v
+
+        PI3D.W = W
+
+        if len(image.shape) == 3:
+            nz,nr,nc = image.shape
+        elif len(image.shape) == 4: # multi-channel image
+            nz,nw,nr,nc = image.shape
+
+        PI3D.NR = nr
+        PI3D.NC = nc
+        PI3D.NZ = nz
+
+        npr = int(np.ceil(nr/subPatchSize)) # number of patch rows
+        npc = int(np.ceil(nc/subPatchSize)) # number of patch cols
+        npz = int(np.ceil(nz/subPatchSize)) # number of patch planes
+
+        nrpi = npr*subPatchSize+2*margin # number of rows in padded image 
+        ncpi = npc*subPatchSize+2*margin # number of cols in padded image 
+        nzpi = npz*subPatchSize+2*margin # number of plns in padded image 
+
+        PI3D.NRPI = nrpi
+        PI3D.NCPI = ncpi
+        PI3D.NZPI = nzpi
+
+        if len(image.shape) == 3:
+            PI3D.PaddedImage = np.zeros((nzpi,nrpi,ncpi))
+            PI3D.PaddedImage[margin:margin+nz,margin:margin+nr,margin:margin+nc] = image
+        elif len(image.shape) == 4:
+            PI3D.PaddedImage = np.zeros((nzpi,nw,nrpi,ncpi))
+            PI3D.PaddedImage[margin:margin+nz,:,margin:margin+nr,margin:margin+nc] = image
+
+        PI3D.PC = [] # patch coordinates [z0,z1,r0,r1,c0,c1]
+        for iZ in range(npz):
+            z0 = iZ*subPatchSize
+            z1 = z0+patchSize
+            for i in range(npr):
+                r0 = i*subPatchSize
+                r1 = r0+patchSize
+                for j in range(npc):
+                    c0 = j*subPatchSize
+                    c1 = c0+patchSize
+                    PI3D.PC.append([z0,z1,r0,r1,c0,c1])
+
+        PI3D.NumPatches = len(PI3D.PC)
+        PI3D.Mode = mode # 'replace' or 'accumulate'
+
+    def getPatch(i):
+        z0,z1,r0,r1,c0,c1 = PI3D.PC[i]
+        if len(PI3D.PaddedImage.shape) == 3:
+            return PI3D.PaddedImage[z0:z1,r0:r1,c0:c1]
+        if len(PI3D.PaddedImage.shape) == 4:
+            return PI3D.PaddedImage[z0:z1,:,r0:r1,c0:c1]
+
+    def createOutput(nChannels):
+        if nChannels == 1:
+            PI3D.Output = np.zeros((PI3D.NZPI,PI3D.NRPI,PI3D.NCPI))
+        else:
+            PI3D.Output = np.zeros((PI3D.NZPI,nChannels,PI3D.NRPI,PI3D.NCPI))
+        if PI3D.Mode == 'accumulate':
+            PI3D.Count = np.zeros((PI3D.NZPI,PI3D.NRPI,PI3D.NCPI))
+
+    def patchOutput(i,P):
+        z0,z1,r0,r1,c0,c1 = PI3D.PC[i]
+        if PI3D.Mode == 'accumulate':
+            PI3D.Count[z0:z1,r0:r1,c0:c1] += PI3D.W
+        if len(P.shape) == 3:
+            if PI3D.Mode == 'accumulate':
+                PI3D.Output[z0:z1,r0:r1,c0:c1] += np.multiply(P,PI3D.W)
+            elif PI3D.Mode == 'replace':
+                PI3D.Output[z0:z1,r0:r1,c0:c1] = P
+        elif len(P.shape) == 4:
+            if PI3D.Mode == 'accumulate':
+                for i in range(P.shape[1]):
+                    PI3D.Output[z0:z1,i,r0:r1,c0:c1] += np.multiply(P[:,i,:,:],PI3D.W)
+            elif PI3D.Mode == 'replace':
+                PI3D.Output[z0:z1,:,r0:r1,c0:c1] = P
+
+    def getValidOutput():
+        margin = PI3D.Margin
+        nz, nr, nc = PI3D.NZ, PI3D.NR, PI3D.NC
+        if PI3D.Mode == 'accumulate':
+            C = PI3D.Count[margin:margin+nz,margin:margin+nr,margin:margin+nc]
+        if len(PI3D.Output.shape) == 3:
+            if PI3D.Mode == 'accumulate':
+                return np.divide(PI3D.Output[margin:margin+nz,margin:margin+nr,margin:margin+nc],C)
+            if PI3D.Mode == 'replace':
+                return PI3D.Output[margin:margin+nz,margin:margin+nr,margin:margin+nc]
+        if len(PI3D.Output.shape) == 4:
+            if PI3D.Mode == 'accumulate':
+                for i in range(PI3D.Output.shape[1]):
+                    PI3D.Output[margin:margin+nz,i,margin:margin+nr,margin:margin+nc] = np.divide(PI3D.Output[margin:margin+nz,i,margin:margin+nr,margin:margin+nc],C)
+            return PI3D.Output[margin:margin+nz,:,margin:margin+nr,margin:margin+nc]
+
+
+    def demo():
+        I = np.random.rand(128,128,128)
+        PI3D.setup(I,64,4,'accumulate')
+
+        nChannels = 2
+        PI3D.createOutput(nChannels)
+
+        for i in range(PI3D.NumPatches):
+            P = PI3D.getPatch(i)
+            Q = np.zeros((P.shape[0],nChannels,P.shape[1],P.shape[2]))
+            for j in range(nChannels):
+                Q[:,j,:,:] = P
+            PI3D.patchOutput(i,Q)
+
+        J = PI3D.getValidOutput()
+        J = J[:,0,:,:]
+
+        D = np.abs(I-J)
+        print(np.max(D))
+
+        pI = I[64,:,:]
+        pJ = J[64,:,:]
+        pD = D[64,:,:]
+
+        K = cat(1,cat(1,pI,pJ),pD)
+        imshow(K)
+