1

I have a 3D image with shape DxHxW. I was successful to extract the image into patches pdxphxpw(overlapping patches). For each patch, I do some processing. Now, I would like to generate the image from the processed patches such that the new image must be same shape with original image. Could you help me to do it.

enter image description here

This is my code to extract patch

def patch_extract_3D(input,patch_shape,xstep=1,ystep=1,zstep=1):
    patches_3D = np.lib.stride_tricks.as_strided(input, ((input.shape[0] - patch_shape[0] + 1) / xstep, (input.shape[1] - patch_shape[1] + 1) / ystep,
                                                  (input.shape[2] - patch_shape[2] + 1) / zstep, patch_shape[0], patch_shape[1], patch_shape[2]),
                                                  (input.strides[0] * xstep, input.strides[1] * ystep,input.strides[2] * zstep, input.strides[0], input.strides[1],input.strides[2]))
    patches_3D= patches_3D.reshape(patches_3D.shape[0]*patches_3D.shape[1]*patches_3D.shape[2], patch_shape[0],patch_shape[1],patch_shape[2])
    return patches_3D

This is the processing the patches (just simple multiple with 2

for i in range(patches_3D.shape[0]):
    patches_3D[i]=patches_3D[i];
    patches_3D[i]=patches_3D[i]*2;

Now, what I need is from patches_3D, I want to reshape it to the original image. Thanks

This is example code

patch_shape=[2, 2, 2]
input=np.arange(4*4*6).reshape(4,4,6)
patches_3D=patch_extract_3D(input,patch_shape)
print  patches_3D.shape
for i in range(patches_3D.shape[0]):
    patches_3D[i]=patches_3D[i]*2
print  patches_3D.shape
user3051460
  • 1,455
  • 22
  • 58
  • `patches_3D` is `pdxphxpw`. An elementwise multiplication of `2` would still keep it as `pdxphxpw`. So, I am not sure how from there you can get to the original image shape of `DxHxW`. – Divakar Feb 24 '17 at 08:41
  • I just multiple intensity, not the shape – user3051460 Feb 24 '17 at 08:43
  • Don't think you understood my point. `patches_3D` and the original image are of different shapes. You can't go back from `patches_3D` to the original one without some sort of reduction. Also, in your sample `patches_3D.shape` is (0,2,2,2). Would make sense if you would test out the samples before posting. – Divakar Feb 24 '17 at 08:55
  • Sorry, the step must be 1 in above test shape 2. the shape of patches_3D is `(45, 2, 2, 2)`. – user3051460 Feb 24 '17 at 09:00

1 Answers1

5

This will do the reverse, however, since your patches overlap this will only be well-defined if their values agree where they overlap

def stuff_patches_3D(out_shape,patches,xstep=12,ystep=12,zstep=12):
    out = np.zeros(out_shape, patches.dtype)
    patch_shape = patches.shape[-3:]
    patches_6D = np.lib.stride_tricks.as_strided(out, ((out.shape[0] - patch_shape[0] + 1) // xstep, (out.shape[1] - patch_shape[1] + 1) // ystep,
                                                  (out.shape[2] - patch_shape[2] + 1) // zstep, patch_shape[0], patch_shape[1], patch_shape[2]),
                                                  (out.strides[0] * xstep, out.strides[1] * ystep,out.strides[2] * zstep, out.strides[0], out.strides[1],out.strides[2]))
    patches_6D[...] = patches.reshape(patches_6D.shape)
    return out

Update: here is a safer version that averages overlapping pixels:

def stuff_patches_3D(out_shape,patches,xstep=12,ystep=12,zstep=12):
    out = np.zeros(out_shape, patches.dtype)
    denom = np.zeros(out_shape, patches.dtype)
    patch_shape = patches.shape[-3:]
    patches_6D = np.lib.stride_tricks.as_strided(out, ((out.shape[0] - patch_shape[0] + 1) // xstep, (out.shape[1] - patch_shape[1] + 1) // ystep,
                                                  (out.shape[2] - patch_shape[2] + 1) // zstep, patch_shape[0], patch_shape[1], patch_shape[2]),
                                                  (out.strides[0] * xstep, out.strides[1] * ystep,out.strides[2] * zstep, out.strides[0], out.strides[1],out.strides[2]))
    denom_6D = np.lib.stride_tricks.as_strided(denom, ((denom.shape[0] - patch_shape[0] + 1) // xstep, (denom.shape[1] - patch_shape[1] + 1) // ystep,
                                                  (denom.shape[2] - patch_shape[2] + 1) // zstep, patch_shape[0], patch_shape[1], patch_shape[2]),
                                                  (denom.strides[0] * xstep, denom.strides[1] * ystep,denom.strides[2] * zstep, denom.strides[0], denom.strides[1],denom.strides[2]))
    np.add.at(patches_6D, tuple(x.ravel() for x in np.indices(patches_6D.shape)), patches.ravel())
    np.add.at(denom_6D, tuple(x.ravel() for x in np.indices(patches_6D.shape)), 1)
    return out/denom
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
  • Thanks for your solution. I just ask about returned variable. Why is out now patches_6D? If returned variable is out, why do you need the `patches_6D[...] = patches.reshape(patches_6D.shape)` – user3051460 Feb 24 '17 at 10:12
  • out and patches_6d reference the same memory, but out has the shape you want. But, please, be really careful with this solution, as it will overwrite the first patches where they overlap with thet later ones. – Paul Panzer Feb 24 '17 at 10:15
  • @user3051460 I've added a safer version. I think it would be better if you used that instead. – Paul Panzer Feb 24 '17 at 10:23
  • Thank you, even it is slower – user3051460 Feb 24 '17 at 10:30
  • Thanks. I would like to ask you about the parameter. Why do we need plus 1 in `(input.shape[0] - patch_shape[0] + 1) / xstep`. I ask this because sometime if we do not setting correct `xstep`, the patch cannot go to near boundary. I guess – user3051460 Feb 24 '17 at 11:44
  • For example, if input is `256x128x256` and patch size is `32x32x32`. Do you think the `xstep,ystep,zstep` can be 16? – user3051460 Feb 24 '17 at 11:45
  • Sure, why not? The logic fiils in as many steps as will fit. Btw. this should be floor divisions (`//` not `/`) I'l change it in the post. For the plus 1 just think of a simple example (patch.shape[0]==1, xstep==1) and you'll see the +1 gives the right number of possible offsets. – Paul Panzer Feb 24 '17 at 12:18