From 8d630fb1b72a57934ab77fc69222b81b2560ea72 Mon Sep 17 00:00:00 2001 From: Tom Runia Date: Wed, 12 Dec 2018 16:15:25 +0100 Subject: [PATCH] debugging --- examples/compare_both.py | 3 + examples/debug_reconstruction.py | 69 +++++++++++++++++ steerable/SCFpyr_NumPy.py | 86 +++++++++++++++------ steerable/SCFpyr_PyTorch.py | 129 +++++++++++++++++-------------- 4 files changed, 205 insertions(+), 82 deletions(-) create mode 100644 examples/debug_reconstruction.py diff --git a/examples/compare_both.py b/examples/compare_both.py index b93f205..7df26cc 100644 --- a/examples/compare_both.py +++ b/examples/compare_both.py @@ -67,6 +67,9 @@ # Also moves the example to CPU and NumPy coeff_torch = utils.extract_from_batch(coeff_torch, 0) +cv2.waitKey(0) +exit() + ################################################################################ # Check correctness diff --git a/examples/debug_reconstruction.py b/examples/debug_reconstruction.py new file mode 100644 index 0000000..ff6e151 --- /dev/null +++ b/examples/debug_reconstruction.py @@ -0,0 +1,69 @@ +# MIT License +# +# Copyright (c) 2018 Tom Runia +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to conditions. +# +# Author: Tom Runia +# Date Created: 2018-12-12 + +from __future__ import absolute_import +from __future__ import division +from __future__ import print_function + +import numpy as np +import cv2 +import cortex.vision + +numpy_outdft = np.load('./assets/numpy_outdft_real.npy') + 1j*np.load('./assets/numpy_outdft_imag.npy') +torch_outdft = np.load('./assets/torch_outdft_real.npy') + 1j*np.load('./assets/torch_outdft_imag.npy') + +print('allclose, real', np.allclose(numpy_outdft.real, torch_outdft.real, atol=1e-2)) +print('allclose, imag', np.allclose(numpy_outdft.imag, torch_outdft.imag, atol=1e-2)) + +real, imag = numpy_outdft.real, numpy_outdft.imag +print(' [numpy] outdft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0:20,0:20].min().item(), real[0:20,0:20].mean().item(), real[0:20,0:20].max().item() +)) +print(' [numpy] outdft imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0:20,0:20].min().item(), imag[0:20,0:20].mean().item(), imag[0:20,0:20].max().item() +)) + +real, imag = torch_outdft.real, torch_outdft.imag +print(' [torch] outdft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0:20,0:20].min().item(), real[0:20,0:20].mean().item(), real[0:20,0:20].max().item() +)) +print(' [torch] outdft imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0:20,0:20].min().item(), imag[0:20,0:20].mean().item(), imag[0:20,0:20].max().item() +)) + +numpy_reconstruction = np.fft.ifftshift(numpy_outdft) +numpy_reconstruction = np.fft.ifft2(numpy_reconstruction) +numpy_reconstruction = numpy_reconstruction.real +real = numpy_reconstruction +print(' [numpy] outdft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0:20,0:20].min().item(), real[0:20,0:20].mean().item(), real[0:20,0:20].max().item() +)) + +torch_reconstruction = np.fft.ifftshift(torch_outdft) +torch_reconstruction = np.fft.ifft2(torch_reconstruction) +torch_reconstruction = torch_reconstruction.real +real = torch_reconstruction +print(' [torch] outdft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0:20,0:20].min().item(), real[0:20,0:20].mean().item(), real[0:20,0:20].max().item() +)) + +cv2.imshow('numpy phase', np.angle(numpy_outdft)) +cv2.imshow('numpy imag', np.abs(numpy_outdft)) +cv2.imshow('torch phase', np.angle(torch_outdft)) +cv2.imshow('torch imag', np.abs(torch_outdft)) + +cv2.imshow('numpy reconstruction', numpy_reconstruction.astype(np.uint8)) +cv2.imshow('torch reconstruction', torch_reconstruction.astype(np.uint8)) + +cv2.waitKey(0) \ No newline at end of file diff --git a/steerable/SCFpyr_NumPy.py b/steerable/SCFpyr_NumPy.py index 52adbfc..a64e769 100644 --- a/steerable/SCFpyr_NumPy.py +++ b/steerable/SCFpyr_NumPy.py @@ -164,20 +164,21 @@ def _build_levels(self, lodft, log_rad, angle, Xrcos, Yrcos, height): return coeff - ################################################################################ - # Reconstruction to Image + ############################################################################ + ########################### RECONSTRUCTION ################################# + ############################################################################ def reconstruct(self, coeff): if self.nbands != len(coeff[1]): raise Exception("Unmatched number of orientations") - M, N = coeff[0].shape - log_rad, angle = math_utils.prepare_grid(M, N) + height, width = coeff[0].shape + log_rad, angle = math_utils.prepare_grid(height, width) Xrcos, Yrcos = math_utils.rcosFn(1, -0.5) Yrcos = np.sqrt(Yrcos) - YIrcos = np.sqrt(np.abs(1 - Yrcos*Yrcos)) + YIrcos = np.sqrt(np.abs(1 - Yrcos**2)) lo0mask = pointOp(log_rad, YIrcos, Xrcos) hi0mask = pointOp(log_rad, Yrcos, Xrcos) @@ -185,25 +186,60 @@ def reconstruct(self, coeff): tempdft = self._reconstruct_levels(coeff[1:], log_rad, Xrcos, Yrcos, angle) hidft = np.fft.fftshift(np.fft.fft2(coeff[0])) + + real, imag = hidft.real, hidft.imag + print(' [numpy] level final. hidft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0:20,0:20].min().item(), real[0:20,0:20].mean().item(), real[0:20,0:20].max().item() + )) + print(' [numpy] level final. hidft imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0:20,0:20].min().item(), imag[0:20,0:20].mean().item(), imag[0:20,0:20].max().item() + )) + outdft = tempdft * lo0mask + hidft * hi0mask - reconstruction = np.fft.ifft2(np.fft.ifftshift(outdft)) - reconstruction = reconstruction.real.astype(int) + real, imag = outdft.real, outdft.imag + np.save('./assets/numpy_outdft_real.npy', real) + np.save('./assets/numpy_outdft_imag.npy', imag) + + print(' [numpy] level final. outdft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0:20,0:20].min().item(), real[0:20,0:20].mean().item(), real[0:20,0:20].max().item() + )) + print(' [numpy] level final. outdft imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0:20,0:20].min().item(), imag[0:20,0:20].mean().item(), imag[0:20,0:20].max().item() + )) + + import cv2 + import cortex.vision + + reconstruction = np.fft.ifftshift(outdft) + + real, imag = reconstruction.real, reconstruction.imag + print(' [numpy] level final. ifftshift real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0:20,0:20].min().item(), real[0:20,0:20].mean().item(), real[0:20,0:20].max().item() + )) + print(' [numpy] level final. ifftshift imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0:20,0:20].min().item(), imag[0:20,0:20].mean().item(), imag[0:20,0:20].max().item() + )) + + real, imag = reconstruction.real, reconstruction.imag + cv2.imshow('numpy real', real) + cv2.imshow('numpy imag', imag) + + reconstruction = np.fft.ifft2(reconstruction) + reconstruction = reconstruction.real + + print(' [numpy] level final. reconstruction real ({:.3f}, {:.3f}, {:.3f})'.format( + reconstruction.min(), reconstruction.mean(), reconstruction.max() + )) + return reconstruction def _reconstruct_levels(self, coeff, log_rad, Xrcos, Yrcos, angle): - - print('[numpy] Call to _reconstruct_levels. remaining = {rem}'.format(rem=len(coeff))) if len(coeff) == 1: - print('[numpy] len(coeff)==1') - print('[numpy] coeff[0].shape', coeff[0].shape, coeff[0].dtype) dft = np.fft.fft2(coeff[0]) - print('[torch] dft after fft', dft.shape, dft.dtype) - tmp = np.fft.fftshift(dft) - print('[torch] dft after fftshift', dft.shape, dft.dtype) - print('[numpy] here!!') - return tmp + dft = np.fft.fftshift(dft) + return dft Xrcos = Xrcos - np.log2(self.scale_factor) @@ -228,6 +264,15 @@ def _reconstruct_levels(self, coeff, log_rad, Xrcos, Yrcos, angle): banddft = np.fft.fftshift(np.fft.fft2(coeff[0][b])) orientdft = orientdft + np.power(np.complex(0, 1), order) * banddft * anglemask * himask + real = orientdft.real + imag = orientdft.imag + print(' [numpy] levels remaining {}. orientdft real ({:.3f}, {:.3f}, {:.3f})'.format( + len(coeff), real.min().item(), real.mean().item(), real.max().item() + )) + print(' [numpy] levels remaining {}. orientdft real ({:.3f}, {:.3f}, {:.3f})'.format( + len(coeff), imag.min().item(), imag.mean().item(), imag.max().item() + )) + #################################################################### ########## Lowpass component are upsampled and convoluted ########## #################################################################### @@ -245,18 +290,9 @@ def _reconstruct_levels(self, coeff, log_rad, Xrcos, Yrcos, angle): nresdft = self._reconstruct_levels(coeff[1:], nlog_rad, Xrcos, Yrcos, nangle) - print(' [numpy] nresdft', nresdft.shape, nresdft.dtype) - print(' [numpy] lomask', lomask.shape, lomask.dtype) - resdft = np.zeros(dims, 'complex') - print(' [numpy] resdft', resdft.shape, resdft.dtype) - print(' [numpy] lostart[0] - loend[0]', lostart[0], loend[0]) - print(' [numpy] lostart[1] - loend[1]', lostart[1], loend[1]) - resdft[lostart[0]:loend[0], lostart[1]:loend[1]] = nresdft * lomask - - return resdft + orientdft diff --git a/steerable/SCFpyr_PyTorch.py b/steerable/SCFpyr_PyTorch.py index ee2d123..f2e0199 100644 --- a/steerable/SCFpyr_PyTorch.py +++ b/steerable/SCFpyr_PyTorch.py @@ -58,7 +58,8 @@ def __init__(self, height=5, nbands=4, scale_factor=2, device=None): self.lutsize = 1024 self.Xcosn = np.pi * np.array(range(-(2*self.lutsize+1), (self.lutsize+2)))/self.lutsize self.alpha = (self.Xcosn + np.pi) % (2*np.pi) - np.pi - self.complex_factor = np.power(np.complex(0, -1), self.nbands - 1) + self.complex_fact_construct = np.power(np.complex(0, -1), self.nbands-1) + self.complex_fact_reconstruct = np.power(np.complex(0, 1), self.nbands-1) ################################################################################ # Construction of Steerable Pyramid @@ -158,8 +159,8 @@ def _build_levels(self, lodft, log_rad, angle, Xrcos, Yrcos, height): # Now multiply with complex number # (x+yi)(u+vi) = (xu-yv) + (xv+yu)i banddft = torch.unbind(banddft, -1) - banddft_real = self.complex_factor.real*banddft[0] - self.complex_factor.imag*banddft[1] - banddft_imag = self.complex_factor.real*banddft[1] + self.complex_factor.imag*banddft[0] + banddft_real = self.complex_fact_construct.real*banddft[0] - self.complex_fact_construct.imag*banddft[1] + banddft_imag = self.complex_fact_construct.real*banddft[1] + self.complex_fact_construct.imag*banddft[0] banddft = torch.stack((banddft_real, banddft_imag), -1) band = math_utils.batch_ifftshift2d(banddft) @@ -202,8 +203,9 @@ def _build_levels(self, lodft, log_rad, angle, Xrcos, Yrcos, height): return coeff - ################################################################################ - # Reconstruction to Image + ############################################################################ + ########################### RECONSTRUCTION ################################# + ############################################################################ def reconstruct(self, coeff): @@ -215,7 +217,7 @@ def reconstruct(self, coeff): Xrcos, Yrcos = math_utils.rcosFn(1, -0.5) Yrcos = np.sqrt(Yrcos) - YIrcos = np.sqrt(np.abs(1 - Yrcos*Yrcos)) + YIrcos = np.sqrt(np.abs(1 - Yrcos**2)) lo0mask = pointOp(log_rad, YIrcos, Xrcos) hi0mask = pointOp(log_rad, Yrcos, Xrcos) @@ -224,19 +226,66 @@ def reconstruct(self, coeff): lo0mask = torch.from_numpy(lo0mask).float()[None,:,:,None].to(self.device) hi0mask = torch.from_numpy(hi0mask).float()[None,:,:,None].to(self.device) - tmp = coeff[1:] + # Start recursive reconstruction tempdft = self._reconstruct_levels(coeff[1:], log_rad, Xrcos, Yrcos, angle) hidft = torch.rfft(coeff[0], signal_ndim=2, onesided=False) hidft = math_utils.batch_fftshift2d(hidft) + + real, imag = torch.unbind(hidft, -1) + print(' [torch] level final. hidft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0,0:20,0:20].min().item(), real[0,0:20,0:20].mean().item(), real[0,0:20,0:20].max().item() + )) + print(' [torch] level final. hidft imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0,0:20,0:20].min().item(), imag[0,0:20,0:20].mean().item(), imag[0,0:20,0:20].max().item() + )) outdft = tempdft * lo0mask + hidft * hi0mask + real, imag = torch.unbind(outdft, -1) + np.save('./assets/torch_outdft_real.npy', real[0].cpu().numpy()) + np.save('./assets/torch_outdft_imag.npy', imag[0].cpu().numpy()) + + print(' [torch] level final. outdft real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0,0:20,0:20].min().item(), real[0,0:20,0:20].mean().item(), real[0,0:20,0:20].max().item() + )) + print(' [torch] level final. outdft imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0,0:20,0:20].min().item(), imag[0,0:20,0:20].mean().item(), imag[0,0:20,0:20].max().item() + )) + # TODO: check these ops (!) - reconstruction = math_utils.batch_fftshift2d(outdft) - reconstruction = torch.ifft(reconstruction, signal_ndim=2) + # hidft = np.fft.fftshift(np.fft.fft2(coeff[0])) + # outdft = tempdft * lo0mask + hidft * hi0mask + + import cv2 + import cortex.vision + + + + # reconstruction = np.fft.ifft2(np.fft.ifftshift(outdft)) + # reconstruction = reconstruction.real.astype(int) + reconstruction = math_utils.batch_ifftshift2d(outdft) + + real, imag = torch.unbind(reconstruction, -1) + print(' [torch] level final. ifftshift real ({:.3f}, {:.3f}, {:.3f})'.format( + real[0,0:20,0:20].min().item(), real[0,0:20,0:20].mean().item(), real[0,0:20,0:20].max().item() + )) + print(' [torch] level final. ifftshift imag ({:.3f}, {:.3f}, {:.3f})'.format( + imag[0,0:20,0:20].min().item(), imag[0,0:20,0:20].mean().item(), imag[0,0:20,0:20].max().item() + )) + + #real, imag = reconstruction, reconstruction #qtorch.unbind(reconstruction, -1) + cv2.imshow('real', real.cpu().numpy()[0]) + cv2.imshow('imag', imag.cpu().numpy()[0]) + + # this fucks it up ... + reconstruction = torch.ifft(reconstruction, signal_ndim=2) #, onesided=False) reconstruction_real = torch.unbind(reconstruction, -1)[0] + print(' [torch] level final. reconstruction real ({:.3f}, {:.3f}, {:.3f})'.format( + reconstruction_real.min().item(), reconstruction_real.mean().item(), reconstruction_real.max().item() + )) + return reconstruction_real def _reconstruct_levels(self, coeff, log_rad, Xrcos, Yrcos, angle): @@ -262,7 +311,6 @@ def _reconstruct_levels(self, coeff, log_rad, Xrcos, Yrcos, angle): Ycosn = np.sqrt(const) * np.power(np.cos(Xcosn), order) orientdft = torch.zeros_like(coeff[0][0]) - for b in range(self.nbands): anglemask = pointOp(angle, Ycosn, Xcosn + np.pi * b/self.nbands) @@ -275,58 +323,34 @@ def _reconstruct_levels(self, coeff, log_rad, Xrcos, Yrcos, angle): # Now multiply with complex number # (x+yi)(u+vi) = (xu-yv) + (xv+yu)i banddft = torch.unbind(banddft, -1) - banddft_real = self.complex_factor.real*banddft[0] - self.complex_factor.imag*banddft[1] - banddft_imag = self.complex_factor.real*banddft[1] + self.complex_factor.imag*banddft[0] + banddft_real = self.complex_fact_reconstruct.real*banddft[0] - self.complex_fact_reconstruct.imag*banddft[1] + banddft_imag = self.complex_fact_reconstruct.real*banddft[1] + self.complex_fact_reconstruct.imag*banddft[0] banddft = torch.stack((banddft_real, banddft_imag), -1) banddft = banddft * anglemask * himask orientdft += banddft + + real, imag = torch.unbind(orientdft, -1) + print(' [torch] levels remaining {}. orientdft real ({:.3f}, {:.3f}, {:.3f})'.format( + len(coeff), real.min().item(), real.mean().item(), real.max().item() + )) + print(' [torch] levels remaining {}. orientdft real ({:.3f}, {:.3f}, {:.3f})'.format( + len(coeff), imag.min().item(), imag.mean().item(), imag.max().item() + )) #################################################################### ########## Lowpass component are upsampled and convoluted ########## #################################################################### - - # Don't consider batch_size and imag/real dim - # dims = np.array(coeff[0][0].shape[1:3]) - - # # Both are tuples of size 2 - # low_ind_start = (np.ceil((dims+0.5)/2) - np.ceil((np.ceil((dims-0.5)/2)+0.5)/2)).astype(int) - # low_ind_end = (low_ind_start + np.ceil((dims-0.5)/2)).astype(int) - - # # Subsampling indices - # nlog_rad = log_rad[low_ind_start[0]:low_ind_end[0],low_ind_start[1]:low_ind_end[1]] - # nangle = angle[low_ind_start[0]:low_ind_end[0],low_ind_start[1]:low_ind_end[1]] - - ################################################################################ - # batch_size = len(coeff[0][0]) - # dims = np.array(coeff[0][0].shape[1:3]) - - # dims = np.array(coeff[0][0].shape[1:3]) - # lostart = (np.ceil((dims+0.5)/2) - np.ceil((np.ceil((dims-0.5)/2)+0.5)/2)).astype(np.int32) - # loend = lostart + np.ceil((dims-0.5)/2).astype(np.int32) - # nlog_rad = log_rad[lostart[0]:loend[0], lostart[1]:loend[1]] # 25,25 - # nangle = angle[lostart[0]:loend[0], lostart[1]:loend[1]] - # YIrcos = np.sqrt(np.abs(1 - Yrcos * Yrcos)) - - # # Filtering - # lomask = pointOp(log_rad, YIrcos, Xrcos) - # lomask = torch.from_numpy(lomask[None,:,:,None]) - # lomask = lomask..float().to(self.device) - - ################################################################################ - # WORKS: - - batch_size = len(coeff[0][0]) dims = np.array(coeff[0][0].shape[1:3]) - print(dims) - lostart = (np.ceil((dims+0.5)/2) - - np.ceil((np.ceil((dims-0.5)/2)+0.5)/2)).astype(np.int32) + lostart = (np.ceil((dims+0.5)/2) - np.ceil((np.ceil((dims-0.5)/2)+0.5)/2)).astype(np.int32) loend = lostart + np.ceil((dims-0.5)/2).astype(np.int32) nlog_rad = log_rad[lostart[0]:loend[0], lostart[1]:loend[1]] nangle = angle[lostart[0]:loend[0], lostart[1]:loend[1]] YIrcos = np.sqrt(np.abs(1 - Yrcos * Yrcos)) + + # Filtering lomask = pointOp(nlog_rad, YIrcos, Xrcos) lomask = torch.from_numpy(lomask[None,:,:,None]) lomask = lomask.float().to(self.device) @@ -337,16 +361,7 @@ def _reconstruct_levels(self, coeff, log_rad, Xrcos, Yrcos, angle): nresdft = self._reconstruct_levels(coeff[1:], nlog_rad, Xrcos, Yrcos, nangle) resdft = torch.zeros_like(coeff[0][0]).to(self.device) - # nresdft is of incorrect size [1,25,13,2] - # lomask of shape [1,25,25,1] - print(' [torch] nresdft', nresdft.shape) - print(' [torch] lomask', lomask.shape) - print(' [torch] resdft', resdft.shape) - print(' [torch] lostart[0] - loend[0]', lostart[0], loend[0]) - print(' [torch] lostart[1] - loend[1]', lostart[1], loend[1]) - - tmp = nresdft * lomask - resdft[:,lostart[0]:loend[0], lostart[1]:loend[1],:] = tmp + resdft[:,lostart[0]:loend[0], lostart[1]:loend[1],:] = nresdft * lomask return resdft + orientdft