# MAD Lab, University at Buffalo # Copyright (C) 2018 Prakhar Jaiswal # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see .



In [ ]:

    
%gui qt

import numpy as np
from numpy import *
import os.path
import subprocess
import binvox_rw
import fftw3f
from mayavi import mlab



In [ ]:

    
def read_voxel(filename):
    fid = open(filename, 'r')
    model = binvox_rw.read_as_3d_array(fid)
    voxel = model.data
    return voxel

def get_sublevel_set(voxel, level):
    return np.nonzero(voxel > 0.9999*level)

def get_volume(voxel):
    coords = get_sublevel_set(voxel, 1)
    return len(coords[0])

def get_fourier_transform(voxel):
    voxel = voxel.astype('f')
    voxel_ft = voxel.astype('F')
    trans = fftw3f.Plan(voxel, voxel_ft, direction = 'forward')
    trans()
    return voxel_ft

def get_inverse_fourier_transform(voxel_ft):
    voxel = zeros(voxel_ft.shape, dtype = 'f')
    trans = fftw3f.Plan(voxel_ft, voxel, direction = 'backward')
    trans()
    return voxel

def get_norm_corr(alpha, beta):
    alpha_ft = get_fourier_transform(alpha)
    beta_ft = get_fourier_transform(beta)
    
    corr_ft = alpha_ft * beta_ft
    
    corr = get_inverse_fourier_transform(corr_ft)
    
    corr /= prod(corr.shape)
    corr = fft.fftshift(corr)
    
    return corr

def minkowski_sum(alpha, beta):
    corr = get_norm_corr(alpha, beta)
    msum = 1 * (corr > 0.0001*(get_volume(beta)-0.5))
    return msum

def minkowski_diff(alpha, beta):
    corr = get_norm_corr(alpha, beta)
    mdiff = 1 * (corr > 1*(get_volume(beta)-0.5))
    return mdiff

def minkowski_sum_and_diff(alpha, beta):
    corr = get_norm_corr(alpha, beta)
    msum = 1 * (corr > 0.0001*(get_volume(beta)-0.5))
    mdiff = 1 * (corr > 1*(get_volume(beta)-0.5))
    return msum, mdiff



In [ ]:

    
def read_rotation_samples(filename):
    samples = np.loadtxt(filename)
    return samples

def get_rot_mat(angles):
    phi, theta, psi = angles
    cphi = cos(phi)
    sphi = sin(phi)
    ctheta = cos(theta)
    stheta = sin(theta)
    cpsi = cos(psi)
    spsi = sin(psi)
    rot = zeros([3, 3])
    rot[0, :] = [cpsi*cphi-spsi*ctheta*sphi, -cpsi*sphi-spsi*ctheta*cphi, spsi*stheta]
    rot[1, :] = [spsi*cphi+cpsi*ctheta*sphi, -spsi*sphi+cpsi*ctheta*cphi, -cpsi*stheta]
    rot[2, :] = [stheta*sphi, stheta*cphi, ctheta]
    return matrix(rot)

def rotate(pts, rot_mat):
    if pts.ndim == 1:
        return squeeze(asarray(dot(rot_mat, pts)))
    return array([squeeze(asarray(dot(rot_mat, point))) for point in pts])

def rotate_voxel(voxel, rot_mat):
    rVoxel = affine_transform(voxel, rot_mat)
    rVoxel =  1 * (rVoxel > 0.1)
    return rVoxel

samples = read_rotation_samples('oim20.eul')
print np.linalg.det(dot(get_rot_mat(samples[1]), get_rot_mat(samples[1])))
print type(get_rot_mat(samples[1]))
pts = array([[1, 2, 3],[1, 2, 3],[1, 2, 3],[1, 2, 3]])
rot = get_rot_mat(samples[233])
print rotate(pts, rot)



In [ ]:

    
input_files = ['data/testpart0.obj', 'data/testpart1.obj']

binvox0 = input_files[0][:input_files[0].rfind('.')] + '.binvox'
binvox1 = input_files[1][:input_files[1].rfind('.')] + '.binvox'

if not os.path.isfile(binvox0):
    subprocess.call('./binvox -d 64 ' + input_files[0], shell = True)
if not os.path.isfile(binvox1):
    subprocess.call('./binvox -d 64 ' + input_files[1], shell = True)



In [ ]:

    
alpha = read_voxel(binvox0)
beta = read_voxel(binvox1)

sza = array(alpha.shape)
szb = array(beta.shape)

szc = sza + szb

dims = [int(pow(2, ceil(log(dim)/log(2)))) for dim in szc]

alpha_pad = zeros(dims, dtype = 'f')
beta_pad = zeros(dims, dtype = 'f')

sid = (dims - sza)/2
eid = sid + sza
alpha_pad[sid[0]:eid[0], sid[1]:eid[1], sid[2]:eid[2]] = alpha

sid = (dims - szb)/2
eid = sid + szb
beta_pad[sid[0]:eid[0], sid[1]:eid[1], sid[2]:eid[2]] = beta

print alpha.shape, alpha_pad.shapes
print beta.shape, beta_pad.shape



In [ ]:

    
msum, mdiff = minkowski_sum_and_diff(alpha_pad, beta_pad)

print get_volume(msum), get_volume(mdiff)

alpha_sf = mlab.pipeline.scalar_field(alpha_pad)
mlab.pipeline.iso_surface(alpha_sf, contours = [1], color = (1, 0, 0), opacity = 0.5)

beta_sf = mlab.pipeline.scalar_field(beta_pad)
mlab.pipeline.iso_surface(beta_sf, contours = [1], color = (0, 0, 1), opacity = 1.0)

msum_sf = mlab.pipeline.scalar_field(msum)
mlab.pipeline.iso_surface(msum_sf, contours = [1], color = (1, 1, 0), opacity = 0.3)

mdiff_sf = mlab.pipeline.scalar_field(mdiff)
mlab.pipeline.iso_surface(mdiff_sf, contours = [1], color = (0, 1, 0), opacity = 0.3)

mlab.show()



In [ ]: