Screen Space Curvature using Cuda and Python

Introduction

This is an example that use Ponca to compute Screen Space Curvature in python using Cuda.

Installation and usage

This example requires the following python packages:

• pycuda
• matplotlib

as well as a working cuda setup (tested with 5.0 and 7.0), and the lastest version of eigen dev branch which manage nvcc compiler.

To compile and run the example, call

cd build && make ssgls_py
cd examples/python && python ssgls.py -s 10 -p ./data/ssgls_sample_wc.png -n ./data/ssgls_sample_normal.png

The script will load the input pictures and compute the curvature for a screenspace neighborhood 10x10 pixels.

Screen-Space Curvature typical input. Left: world coordinates. Right: remapped normal vectors

It will generate this picture

Screen-Space Curvature estimation

as well as the following trace on the standard output

Load CUDA kernel ........... DONE in 0.024 seconds.
Init Input data ............ DONE in 13.390 seconds.
Init queries ............... DONE in 0.720 seconds.
Set memory configuration ... DONE in 0.004 seconds.
Launch kernel .............. DONE in 1.505 seconds.
 
###########Configuration#############
Image size: 1902 x 911
NbQueries:  1732722
Scale:      10
################Results################
Max value:  25.3449363708
Min value:  -31.8370857239
#######################################

Note

Don't worry about the important execution time, this example uses a naive and inefficient approach to represent queries with numpy, making the total process really slow. However, the real computation time in cuda is indicated by the last line of the timing trace. Note that we do not use any re-sampling technique, that is why this step is directly related to the scale parameter.

Cuda programming

Here are the technical details related to the cuda and python biding for screen-space curvature estimation.

Define fitting data structure

class MyPoint
{
public:
    enum {Dim = 3};
    typedef float Scalar;
    typedef Eigen::Matrix<Scalar, Dim, 1> VectorType;
    typedef Eigen::Matrix<Scalar, Dim, Dim> MatrixType;
    typedef Eigen::Matrix<Scalar, 2, 1>   ScreenVectorType;
 
    MULTIARCH inline MyPoint(   const VectorType& _pos        = VectorType::Zero(),
                                const VectorType& _normal     = VectorType::Zero(),
                                const ScreenVectorType& _spos = ScreenVectorType::Zero(),
                                const Scalar _dz = 0.f)
        : m_pos(_pos), m_normal(_normal), m_spos(_spos), m_dz(_dz){}
 
    MULTIARCH inline const VectorType& pos()    const { return m_pos; }
    MULTIARCH inline const VectorType& normal() const { return m_normal; }
    MULTIARCH inline const ScreenVectorType& spos() const { return m_spos; }
    MULTIARCH inline const float & dz() const { return m_dz; }
 
 
    MULTIARCH inline VectorType& pos()   { return m_pos; }
    MULTIARCH inline VectorType& normal()    { return m_normal; }
    MULTIARCH inline ScreenVectorType& spos() { return m_spos; }
    MULTIARCH inline float& dz()     { return m_dz; }
 
 
private:
    ScreenVectorType m_spos;
    VectorType  m_pos, m_normal;
    float m_dz; // depth threshold
};

Define weighting functions

class ProjectWeightFunc: public Grenaille::DistWeightFunc<MyPoint, Grenaille::SmoothWeightKernel<Scalar> >
{
public:
    typedef MyPoint::Scalar Scalar;
    typedef MyPoint::VectorType VectorType;
 
    /*
    Default constructor (needed by Grenaille). Note that the screenspace
    evaluation position is specified as parameter
    */
    MULTIARCH inline ProjectWeightFunc( const Scalar& _t                = 1.f,
                                        const ScreenVectorType& _refPos = ScreenVectorType::Zero(),
                                        const Scalar& _dz               = 0.f)
        : Grenaille::DistWeightFunc<MyPoint, Grenaille::SmoothWeightKernel<Scalar> >(_t), m_refPos(_refPos), m_dz(_dz) {}
 
    MULTIARCH inline Scalar w(const VectorType& _q, const MyPoint&  _attributes) const
    {
        Scalar d  = (_attributes.spos()-m_refPos).norm();
        const Scalar dz = _attributes.dz();
 
        if (d > m_t || dz > m_dz)
            return Scalar(0.);
 
        return m_wk.f(d/m_t);
    }
private:
    ScreenVectorType m_refPos;
    float m_dz;
};

Define fitting primitive

typedef Grenaille::Basket<MyPoint,ProjectWeightFunc,Grenaille::OrientedSphereFit, Grenaille::GLSParam> Gls;

Kernel

extern "C" { 
__global__ void doGLS_kernel(int* _params, //[w, h, scale]
                             float* _positions,
                             float* _normals,
                             float* _result)
{
    
    uint x = (blockIdx.x * blockDim.x) + threadIdx.x;
    uint y = (blockIdx.y * blockDim.y) + threadIdx.y;
    
    const int &width     = _params[0];
    const int &height    = _params[1];
 
    if (x < width && y < height)
    {
        const int &scale = _params[2];
 
        ScreenVectorType refPos;
        refPos << x, y;
 
        int dx, dy; // neighbor offset ids
        int nx, ny; // neighbor ids
 
        Gls gls;
        gls.setWeightFunc(ProjectWeightFunc(scale, refPos));
        gls.init( getVector(x,y,width,height,_positions) );
 
        if (getVector(x,y,width,height,_normals).squaredNorm() == 0.f )
        {
            _result[getId(x,y,width,height,0,1)] = -1.0;
        }
        else
        {
            //_result[getId(x,y,width,height,0,1)] = getVector(x,y,width,height,_normals)(0);
            VectorType p, n;
 
            // collect neighborhood
            VectorType one = VectorType::Zero();
 
            for(dy = -scale; dy != scale; dy++)
            {
                for(dx = -scale; dx != scale; dx++)
                {
                    nx = x+dx;
                    ny = y+dy;
 
 
                    // Check image boundaries
                    if (nx >= 0 && ny >= 0 && nx < width && ny < height)
                    {
                        n = getVector(nx,ny,width,height,_normals);
 
                        // add nei only when the _normal is properly defined
                        // need to use an explicit floating point comparison with pycuda
                        if (n.squaredNorm() != 0.f )
                        {
 
                            // RGB to XYZ remapping
                            n =  2.f * n - one;
                            n.normalize();
 
                            // GLS computation
                            gls.addNeighbor(MyPoint(getVector(nx,ny,width,height,_positions),
                                n,
                                ScreenVectorType(nx,ny)));
                        }
                    }
                }
            }
            // closed form minimization
            gls.finalize();
            _result[getId(x,y,width,height,0,1)] = gls.kappa();
        }
    }
}
 
}

Memory access

We use numpy to format the input data, filled by dimension (in object space) and then by the screen-space coordinates:

__device__ int getId(const int _x,
                     const int _y,
                     const int _width,
                     const int _height,
                     const int _component,
                     const int _nbComponent)
{
    return (_component) + _nbComponent*(_x + _y * _width);
}
 
__device__ VectorType getVector(const int _x,
                                const int _y,
                                const int _width,
                                const int _height,
                                const float* _buffer)
{
    VectorType r;
    r << Scalar(_buffer[getId(_x,_y,_width,_height,0,3)]),
        Scalar(_buffer[getId(_x,_y,_width,_height,1,3)]),
        Scalar(_buffer[getId(_x,_y,_width,_height,2,3)]);
    return r;
}
 
 

Python script

We use numpy to format input data.

# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
import pycuda.autoinit
import pycuda.driver as drv
import numpy
import Image
import time
import sys, getopt
 
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
 
from pycuda.compiler import SourceModule
 
  
  
  
 
start_time=0.0
  
def startTimer( funcname, quiet ):
  global start_time
  
  if not quiet:
    sys.stdout.write(funcname)
    sys.stdout.flush()  
  start_time = time.time()
 
 
def stopTimer(quiet):
  elapsed_time = time.time() - start_time
  if not quiet:
    print 'DONE in {0:.3f} seconds.'.format(elapsed_time)  
 
 
 
 
 
def loadKernel():                   
  module = drv.module_from_file("ssgls.ptx")      
  glsKernel    = module.get_function("doGLS_kernel")
  
  return module, glsKernel
  
  
  
    
 
 
def initInputData(imgNormPath, imgwcPath): 
  i_normals   = Image.open (imgNormPath)
  i_positions = Image.open (imgwcPath)
 
  # Todo: check dimensions for depth, normals and positions
  (w,h)     = i_positions.size     
  length    = w*h;      
 
  normals   = numpy.array(i_normals.getdata()).reshape(3*length,1).astype(numpy.float32) / 255.0;
  positions = numpy.array(i_positions.getdata()).reshape(3*length,1).astype(numpy.float32) / 255.0;
  
  positions = positions * 2.0 - 1.0;
    
  return w, h, normals, positions
  
  
  
    
  
 
def initQueries(w,h):
  nbQueries = w*h
 
  queries   = numpy.zeros(2*nbQueries).astype(numpy.float32);
  for j in range (0,h):
      for i in range (0,w):
          queries[2*(i + j*w)] = i
          queries[2*(i + j*w)+1] = j
  return nbQueries, queries
 
 
 
 
def main(argv):    
    
 
  scale       = 10 # neighborhood size in pixels
  imgNormPath = ''
  imgwcPath   = ''
  quiet       = False
  outputPath  = ''
  
  
  
  try:
    opts, args = getopt.getopt(argv,"s:p:n:o:q",["scale=","positions=","normals=","output=","quiet"])
  except getopt.GetoptError:
    print 'ssgl.py -s <scale> -p <positions> -n <normals>'
    sys.exit(2)
  for opt, arg in opts:
    if opt == '-h':
      print 'ssgl.py -s <scale> -p <positions> -n <normals>'
      sys.exit()
    elif opt in ("-s", "--scale"):
      scale = int(arg)
    elif opt in ("-p", "--imgpath"):
      imgwcPath = arg
    elif opt in ("-n", "--normals"):
      imgNormPath = arg
    elif opt in ("-o", "--output"):
      outputPath = arg
    elif opt in ("-q", "--quiet"):
      quiet = True
      
  
   
  
  
  startTimer("Load CUDA kernel ........... ", quiet)
  module, glsKernel = loadKernel()
  stopTimer(quiet)  
  
  startTimer("Init Input data ............ ", quiet)
  w, h, normals, positions = initInputData (imgNormPath, imgwcPath)
  stopTimer(quiet)    
  
  nbQueries  = w*h
    
  startTimer("Init result array memory ... ", quiet)
  result = numpy.zeros(nbQueries).astype(numpy.float32)
  stopTimer(quiet)
 
 
 
  
  blockSize  = 32
  gridWidth  = w/blockSize
  gridHeight = h/blockSize
  
  w          = numpy.int32(w)
  h          = numpy.int32(h)
  scale      = numpy.int32(scale)
 
  startTimer("Launch kernel .............. ", quiet)
  glsKernel(
          drv.In(numpy.array([w, h, scale])),
          drv.In(positions),
          drv.In(normals),
          drv.Out(result),
          block = (blockSize,  blockSize,  1), 
          grid  = (gridWidth,gridHeight))
          
  stopTimer(quiet)  
  
  
  whereAreNaNs = numpy.isnan(result);
  result[whereAreNaNs] = 0;
 
  if not quiet:
    print '\n###########Configuration#############'
    print 'Image size: {} x {}'.format(w,h)
    print 'NbQueries:  {}'.format(nbQueries)
    print 'Scale:      {}'.format(scale)
    print '################Results################'
    print 'Max value:  {}'.format(result.max())
    print 'Min value:  {}'.format(result.min())
    print '#######################################'
 
  scaleFactor = max(abs(result.max()), abs(result.min())); #/ (2.0*scaleFactor) + 0.5) 
  
  if outputPath == '':
    plt.imshow(result.reshape(h,w), vmin=-scaleFactor, vmax=scaleFactor, cmap='seismic')
    plt.show()
  else:
    plt.imsave(fname=outputPath, arr=result.reshape(h,w), vmin=-1.0, vmax=1.0, cmap='seismic')  
  
if __name__ == "__main__":
  main(sys.argv[1:])