From d7dd95cd9eaa223ff8b86e84e6b1b488ff79bcd5 Mon Sep 17 00:00:00 2001
From: Samuel Fadel <samuelfadel@gmail.com>
Date: Sat, 19 Dec 2015 12:29:33 +0100
Subject: New rendering (VoronoiSplat) -- incomplete.

* Added voronoi-like splatting to points: the same technique from
  Messias et al., (2014)
* It is now possible to change the projection technique during
  runtime (possible, but still requires some work)
---
 skelft.cu | 869 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 869 insertions(+)
 create mode 100644 skelft.cu

(limited to 'skelft.cu')

diff --git a/skelft.cu b/skelft.cu
new file mode 100644
index 0000000..be7662c
--- /dev/null
+++ b/skelft.cu
@@ -0,0 +1,869 @@
+#include <device_functions.h>
+#include "skelft.h"
+#include <cstdio>
+#include <iostream>
+
+
+#include <thrust/host_vector.h> 
+#include <thrust/device_vector.h> 
+#include <thrust/sort.h>
+
+
+
+using namespace std;
+
+
+// Parameters for CUDA kernel executions; more or less optimized for a 1024x1024 image.
+#define BLOCKX     16
+#define BLOCKY     16
+#define BLOCKSIZE  64
+#define TILE_DIM   32
+#define BLOCK_ROWS 16
+
+
+
+/****** Global Variables *******/
+const int NB = 7;                        // Nr buffers we use and store in the entire framework
+short2 **pbaTextures;                    // Work buffers used to compute and store resident images
+// 0: work buffer
+// 1: FT
+// 2: thresholded DT or exact floating-point DT
+// 3: thresholded skeleton
+// 4: topology analysis
+// 5: work buffer for topology
+// 6: skeleton FT
+//
+
+float*  pbaTexSiteParam;    // Stores boundary parameterization
+float2* pbaTexAdvectSites;
+float*  pbaTexDensityData;  // The density map of the sites
+float2* pbaTexGradientData; // The gradient of the density map
+
+int     pbaTexSize;     // Texture size (squared) actually used in all computations
+int     floodBand  = 4, // Various FT computation parameters; defaults are good for an 1024x1024 image.    
+        maurerBand = 4,
+        colorBand  = 4;        
+
+texture<short2> pbaTexColor;  // 2D textures (bound to various buffers defined above as needed)
+texture<short2> pbaTexColor2;
+texture<short2> pbaTexLinks;
+texture<float>  pbaTexParam;  // 1D site parameterization texture (bound to pbaTexSiteParam)
+texture<unsigned char>
+                pbaTexGray;   // 2D texture of unsigned char values, e.g. the binary skeleton
+texture<float2> pbaTexAdvect;
+texture<float,cudaTextureType2D,cudaReadModeElementType>  pbaTexDensity;
+texture<float2,cudaTextureType2D,cudaReadModeElementType> pbaTexGradient;
+texture<float,cudaTextureType2D,cudaReadModeElementType>  pbaTexDT;
+
+
+__device__  bool fill_gc; //Indicates if a fill-sweep did fill anything or not
+
+
+#if __CUDA_ARCH__ < 110   // We cannot use atomic intrinsics on SM10 or below. Thus, we define these as nop.
+#define atomicInc(a,b) 0  // The default will be that some code e.g. endpoint detection will thus not do anything.
+#endif
+
+
+
+/********* Kernels ********/
+#include "skelftkernel.h"
+
+
+
+// Initialize necessary memory (CPU/GPU sides)
+// - textureSize: The max size of any image we will process until re-initialization
+void skelft2DInitialization(int maxTexSize)
+{
+    int pbaMemSize = maxTexSize * maxTexSize * sizeof(short2); // A buffer has 2 shorts / pixel
+    pbaTextures  = (short2**) malloc(NB * sizeof(short2*));    // We will use NB buffers
+
+    for(int i=0;i<NB;++i)
+       cudaMalloc((void **) &pbaTextures[i], pbaMemSize); // Allocate work buffer 'i' (on GPU)
+
+    cudaMalloc((void**) &pbaTexSiteParam, maxTexSize * maxTexSize * sizeof(float));    // Sites texture (for FT)
+    cudaMalloc((void**) &pbaTexAdvectSites, maxTexSize * maxTexSize * sizeof(float2)); // Sites 2D-coords
+    cudaMalloc((void**) &pbaTexDensityData, maxTexSize * maxTexSize * sizeof(float));
+    cudaMalloc((void**) &pbaTexGradientData, maxTexSize * maxTexSize * sizeof(float2));
+}
+
+
+
+
+// Deallocate all allocated memory
+void skelft2DDeinitialization()
+{
+    for(int i=0;i<NB;++i) cudaFree(pbaTextures[i]); 
+    cudaFree(pbaTexSiteParam);
+    cudaFree(pbaTexAdvectSites);
+    cudaFree(pbaTexDensityData);
+    cudaFree(pbaTexGradientData);
+    free(pbaTextures);
+}
+
+
+// Initialize the Voronoi textures from the sites' encoding texture (parameterization)
+// REMARK: we interpret 'inputVoro' as a 2D texture, as it's much easier/faster like this
+__global__ void kernelSiteParamInit(short2* inputVoro, int size)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+
+    if (tx<size && ty<size) // Careful not to go outside the image..
+    {
+        int i = TOID(tx,ty,size);
+        // The sites-param has non-zero (parameter) values precisely on boundary points
+        float param = tex1Dfetch(pbaTexParam,i); 
+
+        short2& v = inputVoro[i];
+        // Non-boundary points are marked as 0 in the parameterization. Here we will compute the FT.
+        v.x = v.y = MARKER;
+        // These are points which define the 'sites' to compute the FT/skeleton (thus, have FT==identity)
+        // We could use an if-then-else here, but it's faster with an if-then
+        if (param>0) 
+        {            
+            v.x = tx;
+            v.y = ty;
+        }
+    }
+}
+
+
+
+void skelft2DInitializeInput(float* sites, int size)                                        // Copy input sites from CPU to GPU; Also set up site param initialization in pbaTextures[0]
+{
+    pbaTexSize = size;                                                                        // Size of the actual texture being used in this run; can be smaller than the max-tex-size
+                                                                                            // which was used in skelft2DInitialization()
+
+    cudaMemcpy(pbaTexSiteParam, sites, pbaTexSize * pbaTexSize * sizeof(float), cudaMemcpyHostToDevice);
+                                                                                            // Pass sites parameterization to CUDA.  Must be done before calling the initialization
+                                                                                            // kernel, since we use the sites-param as a texture in that kernel
+    cudaBindTexture(0, pbaTexParam, pbaTexSiteParam);                                        // Bind the sites-param as a 1D texture so we can quickly index it next
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+    
+    kernelSiteParamInit<<<grid,block>>>(pbaTextures[0],pbaTexSize);                            // Do the site param initialization. This sets up pbaTextures[0]
+    cudaUnbindTexture(pbaTexParam);
+}
+
+
+
+
+
+// In-place transpose a squared texture. 
+// Block orders are modified to optimize memory access. 
+// Point coordinates are also swapped. 
+void pba2DTranspose(short2 *texture)
+{
+    dim3 block(TILE_DIM, BLOCK_ROWS); 
+    dim3 grid(pbaTexSize / TILE_DIM, pbaTexSize / TILE_DIM); 
+
+    cudaBindTexture(0, pbaTexColor, texture); 
+    kernelTranspose<<< grid, block >>>(texture, pbaTexSize); 
+    cudaUnbindTexture(pbaTexColor); 
+}
+
+// Phase 1 of PBA. m1 must divides texture size
+void pba2DPhase1(int m1, short xm, short ym, short xM, short yM) 
+{
+    dim3 block = dim3(BLOCKSIZE);   
+    dim3 grid = dim3(pbaTexSize / block.x, m1); 
+
+    // Flood vertically in their own bands
+    cudaBindTexture(0, pbaTexColor, pbaTextures[0]); 
+    kernelFloodDown<<< grid, block >>>(pbaTextures[1], pbaTexSize, pbaTexSize / m1); 
+    cudaUnbindTexture(pbaTexColor); 
+
+    cudaBindTexture(0, pbaTexColor, pbaTextures[1]); 
+    kernelFloodUp<<< grid, block >>>(pbaTextures[1], pbaTexSize, pbaTexSize / m1); 
+
+    // Passing information between bands
+    grid = dim3(pbaTexSize / block.x, m1); 
+    kernelPropagateInterband<<< grid, block >>>(pbaTextures[0], pbaTexSize, pbaTexSize / m1); 
+
+    cudaBindTexture(0, pbaTexLinks, pbaTextures[0]); 
+    kernelUpdateVertical<<< grid, block >>>(pbaTextures[1], pbaTexSize, m1, pbaTexSize / m1); 
+    cudaUnbindTexture(pbaTexLinks); 
+    cudaUnbindTexture(pbaTexColor); 
+}
+
+// Phase 2 of PBA. m2 must divides texture size
+void pba2DPhase2(int m2) 
+{
+    // Compute proximate points locally in each band
+    dim3 block = dim3(BLOCKSIZE);   
+    dim3 grid = dim3(pbaTexSize / block.x, m2); 
+    cudaBindTexture(0, pbaTexColor, pbaTextures[1]); 
+    kernelProximatePoints<<< grid, block >>>(pbaTextures[0], pbaTexSize, pbaTexSize / m2); 
+
+    cudaBindTexture(0, pbaTexLinks, pbaTextures[0]); 
+    kernelCreateForwardPointers<<< grid, block >>>(pbaTextures[0], pbaTexSize, pbaTexSize / m2); 
+
+    // Repeatly merging two bands into one
+    for (int noBand = m2; noBand > 1; noBand /= 2) {
+        grid = dim3(pbaTexSize / block.x, noBand / 2); 
+        kernelMergeBands<<< grid, block >>>(pbaTextures[0], pbaTexSize, pbaTexSize / noBand); 
+    }
+
+    // Replace the forward link with the X coordinate of the seed to remove
+    // the need of looking at the other texture. We need it for coloring.
+    grid = dim3(pbaTexSize / block.x, pbaTexSize); 
+    kernelDoubleToSingleList<<< grid, block >>>(pbaTextures[0], pbaTexSize); 
+    cudaUnbindTexture(pbaTexLinks); 
+    cudaUnbindTexture(pbaTexColor); 
+}
+
+// Phase 3 of PBA. m3 must divides texture size
+void pba2DPhase3(int m3) 
+{
+    dim3 block = dim3(BLOCKSIZE / m3, m3); 
+    dim3 grid = dim3(pbaTexSize / block.x); 
+    cudaBindTexture(0, pbaTexColor, pbaTextures[0]); 
+    kernelColor<<< grid, block >>>(pbaTextures[1], pbaTexSize); 
+    cudaUnbindTexture(pbaTexColor); 
+}
+
+
+
+void skel2DFTCompute(short xm, short ym, short xM, short yM, int floodBand, int maurerBand, int colorBand)
+{
+    pba2DPhase1(floodBand,xm,ym,xM,yM);                                        //Vertical sweep
+
+    pba2DTranspose(pbaTextures[1]);                                            //
+
+    pba2DPhase2(maurerBand);                                                //Horizontal coloring
+
+    pba2DPhase3(colorBand);                                                    //Row coloring
+
+    pba2DTranspose(pbaTextures[1]); 
+}
+
+
+
+
+
+__global__ void kernelThresholdDT(unsigned char* output, int size, float threshold2, short xm, short ym, short xM, short yM)
+//Input:    pbaTexColor: closest-site-ids per pixel, i.e. FT
+//Output:   output: thresholded DT
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+    
+    if (tx>xm && ty>ym && tx<xM && ty<yM)                                    //careful not to index outside the image..
+    {    
+        int    id     = TOID(tx, ty, size);
+      short2 voroid = tex1Dfetch(pbaTexColor,id);                            //get the closest-site to tx,ty into voroid.x,.y
+      float  d2     = (tx-voroid.x)*(tx-voroid.x)+(ty-voroid.y)*(ty-voroid.y);
+      output[id]    = (d2<=threshold2);                                        //threshold DT into binary image
+    }
+}
+
+
+
+__global__ void kernelDT(float* output, int size, short xm, short ym, short xM, short yM)
+//Input:    pbaTexColor: closest-site-ids per pixel, i.e. FT
+//Output:   output: DT
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+        
+    if (tx>xm && ty>ym && tx<xM && ty<yM)                                    //careful not to index outside the image..
+    {    
+        int    id     = TOID(tx, ty, size);
+      short2 voroid = tex1Dfetch(pbaTexColor,id);                            //get the closest-site to tx,ty into voroid.x,.y
+      float  d2     = (tx-voroid.x)*(tx-voroid.x)+(ty-voroid.y)*(ty-voroid.y);
+      output[id]    = sqrtf(d2);                                            //save the Euclidean DT
+    }
+}
+
+
+__global__ void kernelSkel(unsigned char* output, short xm, short ym, 
+                           short xM, short yM, short size, float threshold, float length)    
+                                                                            //Input:    pbaTexColor: closest-site-ids per pixel
+                                                                            //            pbaTexParam: labels for sites (only valid at site locations)
+{                                                                            //Output:    output: binary thresholded skeleton
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+    
+    if (tx>xm && ty>ym && tx<xM && ty<yM) 
+    {
+        int    id     = TOID(tx, ty, size);
+      int    Id     = id;
+      short2 voroid = tex1Dfetch(pbaTexColor,id);                            //get the closest-site to tx,ty into voroid.x,.y
+      int    id2    = TOID(voroid.x,voroid.y,size);                            //convert the site's coord to an index into pbaTexParam[], the site-label-texture
+      float  imp    = tex1Dfetch(pbaTexParam,id2);                            //get the site's label
+
+             ++id;                                                            //TOID(tx+1,ty,size)
+             voroid = tex1Dfetch(pbaTexColor,id);                            //
+             id2    = TOID(voroid.x,voroid.y,size);                            //
+      float  imp_r  = tex1Dfetch(pbaTexParam,id2);                            //
+
+             id     += size-1;                                                //TOID(tx,ty+1,size)
+             voroid = tex1Dfetch(pbaTexColor,id);                            //
+             id2    = TOID(voroid.x,voroid.y,size);                            //
+      float  imp_u  = tex1Dfetch(pbaTexParam,id2);                            //
+
+      float imp_dx  = fabsf(imp_r-imp);
+      float imp_dy  = fabsf(imp_u-imp);
+      float Imp     = max(imp_dx,imp_dy);
+      Imp = min(Imp,fabsf(length-Imp));
+      if (Imp>=threshold) output[Id] = 1;                                    //By filling only 1-values, we reduce memory access somehow (writing to output[] is expensive)
+    } 
+    
+    //WARNING: this kernel may sometimes creates 2-pixel-thick branches.. Study the AFMM original code to see if this is correct.
+}
+
+
+
+#define X 1
+
+__constant__ const                                                            //REMARK: put following constants (for kernelTopology) in CUDA constant-memory, as this gives a huge speed difference
+unsigned char topo_patterns[][9] =        { {0,0,0,                            //These are the 3x3 templates that we use to detect skeleton endpoints
+                                           0,X,0,                            //(with four 90-degree rotations for each)
+                                           0,X,0},
+                                          {0,0,0,
+                                           0,X,0,
+                                           0,0,X},
+                                          {0,0,0,
+                                           0,X,0,
+                                           0,X,X},
+                                          {0,0,0,
+                                           0,X,0,
+                                           X,X,0} 
+                                        };
+                                        
+#define topo_NPATTERNS  4                                                        //Number of patterns we try to match (for kernelTopology)
+                                                                                //REMARK: #define faster than __constant__
+
+__constant__ const unsigned char topo_rot[][9] = { {0,1,2,3,4,5,6,7,8}, {2,5,8,1,4,7,0,3,6}, {8,7,6,5,4,3,2,1,0}, {6,3,0,7,4,1,8,5,2} };
+                                                                                //These encode the four 90-degree rotations of the patterns (for kernelTopology);
+
+__device__ unsigned int topo_gc            = 0;
+__device__ unsigned int topo_gc_last    = 0;
+
+
+__global__ void kernelTopology(unsigned char* output, short2* output_set, short xm, short ym, short xM, short yM, short size, int maxpts)    
+{
+    const int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    const int ty = blockIdx.y * blockDim.y + threadIdx.y;
+        
+    unsigned char t[9];
+    
+    if (tx>xm && ty>ym && tx<xM-1 && ty<yM-1)                                    //careful not to index outside the image; take into account the template size too
+    {    
+       int    id = TOID(tx, ty, size);     
+       unsigned char  p  = tex1Dfetch(pbaTexGray,id);                            //get the skeleton pixel at tx,ty
+       if (p)                                                                    //if the pixel isn't skeleton, nothing to do
+       {
+         unsigned char idx=0;
+         for(int j=ty-1;j<=ty+1;++j)                                            //read the template into t[] for easier use
+         {
+           int id = TOID(tx-1, j, size);
+           for(int i=0;i<=2;++i,++id,++idx)
+              t[idx] = tex1Dfetch(pbaTexGray,id);                                //get the 3x3 template centered at the skel point tx,ty
+         }
+          
+         for(unsigned char r=0;r<4;++r)                                            //try to match all rotations of a pattern:
+         {
+           const unsigned char* rr = topo_rot[r];
+           for(unsigned char p=0;p<topo_NPATTERNS;++p)                            //try to match all patterns:
+           {
+             const unsigned char* pat = topo_patterns[p];
+             unsigned char j = (p==0)? 0 : 7;                                    //Speedup: for all patterns except 1st, check only last 3 entries, the first 6 are identical for all patterns
+             for(;j<9;++j)                                                        //try to match rotated pattern vs actual pattern
+               if (pat[j]!=t[rr[j]]) break;                                        //this rotation failed
+             if (j<6) break;                                                    //Speedup: if we have a mismatch on the 1st 6 pattern entries, then none of the patterns can match
+                                                                                //           since all templates have the same first 6 entries.
+
+             if (j==9)                                                            //this rotation succeeded: mark the pixel as a topology event and we're done
+             {    
+                int crt_gc = atomicInc(&topo_gc,maxpts);                        //REMARK: this serializes (compacts) all detected endpoints in one array.            
+                output_set[crt_gc] = make_short2(tx,ty);                        //To do this, we use an atomic read-increment-return on a global counter, 
+                                                                                //which is guaranteed to give all threads unique consecutive indexes in the array.
+                output[id] = 1;                                                    //Also create the topology image
+                return;
+             }
+           }
+         }
+       }
+    }
+    else                                                                        //Last thread: add zero-marker to the output point-set, so the reader knows how many points are really in there
+    if (tx==xM-1 && ty==yM-1)                                                    //Also reset the global vector counter topo_gc, for the next parallel-run of this function
+    { topo_gc_last = topo_gc; topo_gc = 0; }                                    //We do this in the last thread so that no one modifies topo_gc from now on. 
+                                                                                //REMARK: this seems to be the only way I can read a __device__ variable back to the CPU
+}
+
+
+
+
+void skelft2DParams(int floodBand_, int maurerBand_, int colorBand_)        //Set up some params of the FT algorithm    
+{
+  floodBand   = floodBand_;
+  maurerBand  = maurerBand_;
+  colorBand   = colorBand_;
+}
+
+
+
+
+
+// Compute 2D FT / Voronoi diagram of a set of sites
+// siteParam:   Site parameterization. 0 = non-site points; >0 = site parameter value.
+// output:        FT. The (x,y) at (i,j) are the coords of the closest site to (i,j)
+// size:        Texture size (pow 2)
+void skelft2DFT(short* output, float* siteParam, short xm, short ym, short xM, short yM, int size) 
+{
+    skelft2DInitializeInput(siteParam,size);                                    // Initialization of already-allocated data structures
+
+    skel2DFTCompute(xm, ym, xM, yM, floodBand, maurerBand, colorBand);            // Compute FT
+    
+    // Copy FT to CPU, if required
+    if (output) cudaMemcpy(output, pbaTextures[1], size*size*sizeof(short2), cudaMemcpyDeviceToHost);
+}
+
+
+
+
+
+
+void skelft2DDT(float* outputDT, short xm, short ym, short xM, short yM)        //Compute exact DT from resident FT (in pbaTextures[1])
+{
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+    
+
+    cudaBindTexture(0, pbaTexColor, pbaTextures[1]);                            //Used to read the FT from
+
+    xm = ym = 0; xM = yM = pbaTexSize-1;
+    kernelDT <<< grid, block >>>((float*)pbaTextures[2], pbaTexSize, xm-1, ym-1, xM+1, yM+1);
+    cudaUnbindTexture(pbaTexColor);
+      
+    //Copy DT to CPU
+    if (outputDT) cudaMemcpy(outputDT, pbaTextures[2], pbaTexSize * pbaTexSize * sizeof(float), cudaMemcpyDeviceToHost);
+}
+
+
+
+void skelft2DDT(short* outputDT, float threshold,                                //Compute (thresholded) DT (into pbaTextures[2]) from resident FT (in pbaTextures[1])    
+                      short xm, short ym, short xM, short yM)
+{
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+
+    cudaBindTexture(0, pbaTexColor, pbaTextures[1]);                            //Used to read the FT from
+
+    if (threshold>=0)
+    {
+      xm -= (short)threshold; if (xm<0) xm=0;
+      ym -= (short)threshold; if (ym<0) ym=0;
+      xM += (short)threshold; if (xM>pbaTexSize-1) xM=pbaTexSize-1;
+      yM += (short)threshold; if (yM>pbaTexSize-1) yM=pbaTexSize-1;
+    
+      kernelThresholdDT<<< grid, block >>>((unsigned char*)pbaTextures[2], pbaTexSize, threshold*threshold, xm-1, ym-1, xM+1, yM+1);    
+      cudaUnbindTexture(pbaTexColor);
+    
+      //Copy thresholded image to CPU
+      if (outputDT) cudaMemcpy(outputDT, (unsigned char*)pbaTextures[2], pbaTexSize * pbaTexSize * sizeof(unsigned char), cudaMemcpyDeviceToHost);
+    }
+}
+
+
+
+
+void skelft2DSkeleton(unsigned char* outputSkel, float length, float threshold,    //Compute thresholded skeleton (into pbaTextures[3]) from resident FT (in pbaTextures[1])
+                      short xm,short ym,short xM,short yM)                        
+{                                                                                //length:     boundary length
+    dim3 block = dim3(BLOCKX,BLOCKY);                                            //threshold:  skeleton importance min-value (below this, we ignore branches)
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+    
+    cudaBindTexture(0, pbaTexColor, pbaTextures[1]);                            //Used to read the resident FT
+    cudaBindTexture(0, pbaTexParam, pbaTexSiteParam);                            //Used to read the resident boundary parameterization    
+    cudaMemset(pbaTextures[3],0,sizeof(unsigned char)*pbaTexSize*pbaTexSize);    //Faster to zero result and then fill only 1-values (see kernel)
+    
+    kernelSkel<<< grid, block >>>((unsigned char*)pbaTextures[3], xm, ym, xM-1, yM-1, pbaTexSize, threshold, length);
+    
+    cudaUnbindTexture(pbaTexColor);
+    cudaUnbindTexture(pbaTexParam);
+    
+    //Copy skeleton to CPU
+    if (outputSkel) cudaMemcpy(outputSkel, pbaTextures[3], pbaTexSize * pbaTexSize * sizeof(unsigned char), cudaMemcpyDeviceToHost);
+}
+
+
+
+
+void skelft2DTopology(unsigned char* outputTopo, int* npts, short* outputPoints, //Compute topology-points of the resident skeleton (in pbaTextures[3])
+                      short xm,short ym,short xM,short yM)                    
+{
+    int maxpts = (npts)? *npts : pbaTexSize*pbaTexSize;                            //This is the max # topo-points we are going to return in outputPoints[]
+
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+    
+    cudaBindTexture(0, pbaTexGray, pbaTextures[3]);                                //Used to read the resident skeleton
+    cudaMemset(pbaTextures[4],0,sizeof(unsigned char)*pbaTexSize*pbaTexSize);    //Faster to zero result and then fill only 1-values (see kernel)
+
+    unsigned int zero = 0;
+    cudaMemcpyToSymbol(topo_gc,&zero,sizeof(unsigned int),0,cudaMemcpyHostToDevice);        //Set topo_gc to 0
+
+    kernelTopology<<< grid, block >>>((unsigned char*)pbaTextures[4], pbaTextures[5], xm, ym, xM, yM, pbaTexSize, maxpts+1);
+    
+    cudaUnbindTexture(pbaTexGray);
+
+    if (outputPoints && maxpts)                                                    //If output-point vector desired, copy the end-points, put in pbaTexture[5] as a vector of short2's, 
+    {                                                                            //into caller space. We copy only 'maxpts' elements, as the user instructed us.
+        unsigned int num_pts;
+        cudaMemcpyFromSymbol(&num_pts,topo_gc_last,sizeof(unsigned int),0,cudaMemcpyDeviceToHost);        //Get #topo-points we have detected from the device-var from CUDA
+        if (npts && num_pts)                                                                            //Copy the topo-points to caller        
+           cudaMemcpy(outputPoints,pbaTextures[5],num_pts*sizeof(short2),cudaMemcpyDeviceToHost);
+        if (npts) *npts = num_pts;                                                //Return #detected topo-points to caller                                
+    }
+        
+    if (outputTopo)                                                                //If topology image desired, copy it into user space
+        cudaMemcpy(outputTopo,pbaTextures[4],pbaTexSize*pbaTexSize*sizeof(unsigned char), cudaMemcpyDeviceToHost);
+}
+
+
+
+
+__global__ void kernelSiteFromSkeleton(short2* outputSites, int size)                        //Initialize the Voronoi textures from the sites' encoding texture (parameterization)
+{                                                                                            //REMARK: we interpret 'inputVoro' as a 2D texture, as it's much easier/faster like this
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+
+    if (tx<size && ty<size)                                                                    //Careful not to go outside the image..
+    {
+      int i = TOID(tx,ty,size);
+      unsigned char param = tex1Dfetch(pbaTexGray,i);                                        //The sites-param has non-zero (parameter) values precisely on non-boundary points
+
+      short2& v = outputSites[i];
+      v.x = v.y = MARKER;                                                                    //Non-boundary points are marked as 0 in the parameterization. Here we will compute the FT.
+      if (param)                                                                            //These are points which define the 'sites' to compute the FT/skeleton (thus, have FT==identity)
+      {                                                                                        //We could use an if-then-else here, but it's faster with an if-then
+         v.x = tx; v.y = ty;
+      }
+    }
+}
+
+
+
+
+__global__ void kernelSkelInterpolate(float* output, int size)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+
+    if (tx<size && ty<size)                                                                    //Careful not to go outside the image..
+    {
+        int    id     = TOID(tx, ty, size);
+      short2 vid    = tex1Dfetch(pbaTexColor,id);                            
+      float  T      = sqrtf((tx-vid.x)*(tx-vid.x)+(ty-vid.y)*(ty-vid.y));
+      short2 vid2   = tex1Dfetch(pbaTexColor2,id);                            
+      float  D      = sqrtf((tx-vid2.x)*(tx-vid2.x)+(ty-vid2.y)*(ty-vid2.y));
+      float  B      = ((D)? min(T/2/D,0.5f):0.5) + 0.5*((T)? max(1-D/T,0.0f):0);
+      output[id]    = pow(B,0.2f);
+    }
+}
+
+
+
+
+void skel2DSkeletonDT(float* outputSkelDT,short xm,short ym,short xM,short yM)
+{
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+
+    cudaBindTexture(0,pbaTexGray,pbaTextures[3]);                            //Used to read the resident binary skeleton
+    kernelSiteFromSkeleton<<<grid,block>>>(pbaTextures[0],pbaTexSize);        //1. Init pbaTextures[0] with sites on skeleton i.e. from pbaTexGray
+    cudaUnbindTexture(pbaTexGray);
+        
+    //Must first save pbaTextures[1] since we may need it later..
+    cudaMemcpy(pbaTextures[5],pbaTextures[1],pbaTexSize*pbaTexSize*sizeof(short2),cudaMemcpyDeviceToDevice);
+    skel2DFTCompute(xm, ym, xM, yM, floodBand, maurerBand, colorBand);        //2. Compute FT of the skeleton into pbaTextures[6]
+    cudaMemcpy(pbaTextures[6],pbaTextures[1],pbaTexSize*pbaTexSize*sizeof(short2),cudaMemcpyDeviceToDevice);
+    cudaMemcpy(pbaTextures[1],pbaTextures[5],pbaTexSize*pbaTexSize*sizeof(short2),cudaMemcpyDeviceToDevice);
+    
+    //Compute interpolation        
+    cudaBindTexture(0,pbaTexColor,pbaTextures[1]);                            // FT of boundary
+    cudaBindTexture(0,pbaTexColor2,pbaTextures[6]);                            // FT of skeleton
+    kernelSkelInterpolate<<<grid,block>>>((float*)pbaTextures[0],pbaTexSize);
+    cudaUnbindTexture(pbaTexColor);
+    cudaUnbindTexture(pbaTexColor2);
+    if (outputSkelDT) cudaMemcpy(outputSkelDT, pbaTextures[0], pbaTexSize * pbaTexSize * sizeof(float), cudaMemcpyDeviceToHost);
+}
+
+
+
+
+
+__global__ void kernelFill(unsigned char* output, int size, unsigned char bg, unsigned char fg, short xm, short ym, short xM, short yM, bool ne)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+        
+    if (tx>xm && ty>ym && tx<xM && ty<yM)                                        //careful not to index outside the image..
+    {    
+      int    id0 = TOID(tx, ty, size);
+      unsigned char val = tex1Dfetch(pbaTexGray,id0);                            //
+      if (val==fg)                                                                //do we have a filled pixel? Then fill all to left/top/up/bottom of it which is background
+      {
+        bool fill = false;
+        int id = id0;
+        if (ne)                                                                    //fill in north-east direction:
+        {
+            for(short x=tx+1;x<xM;++x)                                            //REMARK: here and below, the interesting thing is that it's faster, by about 10-15%, to fill a whole
+            {                                                                    //        scanline rather than oly until the current block's borders (+1). The reason is that filling a whole
+                                                                                //          scanline decreases the total #sweeps, which seems to be the limiting speed factor
+              if (tex1Dfetch(pbaTexGray,++id)!=bg) break;
+              output[id] = fg; fill = true;
+            }
+
+            id = id0;
+            for(short y=ty-1;y>ym;--y)
+            {
+              if (tex1Dfetch(pbaTexGray,id-=size)!=bg) break;
+              output[id] = fg; fill = true;
+            }
+        }
+        else                                                                    //fill in south-west direction:
+        {
+            for(short x=tx-1;x>xm;--x)
+            {
+              if (tex1Dfetch(pbaTexGray,--id)!=bg) break;
+              output[id] = fg; fill = true;
+            }
+
+            id = id0;
+            for(short y=ty+1;y<yM;++y)
+            {
+              if (tex1Dfetch(pbaTexGray,id+=size)!=bg) break;
+              output[id] = fg; fill = true;
+            }
+        }
+        
+        if (fill) fill_gc = true;                                                //if we filled anything, inform caller; we 'gather' this info from a local var into the
+                                                                                //global var here, since it's faster than writing the global var in the for loops
+      }      
+    }
+}
+
+
+
+
+__global__ void kernelFillHoles(unsigned char* output, int size, unsigned char bg, unsigned char fg, unsigned char fill_fg)
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+    
+    if (tx>=0 && ty>=0 && tx<size && ty<size)                                    //careful not to index outside the image..
+    {    
+        int            id = TOID(tx, ty, size);
+      unsigned char val = tex1Dfetch(pbaTexGray,id);                            //
+      if (val==fill_fg)
+         output[id] = bg;
+      else if (val==bg)
+         output[id] = fg;     
+    }
+}
+
+
+int skelft2DFill(unsigned char* outputFill, short sx, short sy, short xm, short ym, short xM, short yM, unsigned char fill_value)
+{
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+
+    unsigned char background;
+    int id = sy * pbaTexSize + sx;
+    cudaMemcpy(&background,(unsigned char*)pbaTextures[2]+id,sizeof(unsigned char),cudaMemcpyDeviceToHost); //See which is the value we have to fill from (sx,sy)
+    
+    cudaMemset(((unsigned char*)pbaTextures[2])+id,fill_value,sizeof(unsigned char));                    //Fill the seed (x,y) on the GPU    
+
+    cudaBindTexture(0, pbaTexGray, pbaTextures[2]);                                                        //Used to read the thresholded DT
+
+    int iter=0;
+    bool xy = true;                                                                                        //Direction of filling for current sweep: either north-east or south-west
+                                                                                                        //This kind of balances the memory-accesses nicely over kernel calls
+    for(;;++iter,xy=!xy)                                                                                //Keep filling a sweep at a time until we have no background pixels anymore
+    {    
+       bool filled = false;                                                                                //Initialize flag: we didn't fill anything in this sweep
+       cudaMemcpyToSymbol(fill_gc,&filled,sizeof(bool),0,cudaMemcpyHostToDevice);                        //Pass flag to CUDA
+       kernelFill<<<grid, block>>>((unsigned char*)pbaTextures[2],pbaTexSize,background,fill_value,xm,ym,xM,yM,xy);    
+                                                                                                        //One fill sweep       
+       cudaMemcpyFromSymbol(&filled,fill_gc,sizeof(bool),0,cudaMemcpyDeviceToHost);                        //See if we filled anything in this sweep
+       if (!filled) break;                                                                                //Nothing filled? Then we're done, the image didn't change
+    }
+    cudaUnbindTexture(pbaTexGray);
+        
+    if (outputFill) cudaMemcpy(outputFill, (unsigned char*)pbaTextures[2], pbaTexSize * pbaTexSize * sizeof(unsigned char), cudaMemcpyDeviceToHost);
+    
+    return iter;                                                                                        //Return #iterations done for the fill - useful as a performance measure for caller
+}
+
+
+
+int skelft2DFillHoles(unsigned char* outputFill, short sx, short sy, unsigned char foreground)
+{
+    unsigned char background;
+    unsigned char fill_value = 128;
+    int id = sy * pbaTexSize + sx;
+    cudaMemcpy(&background,(unsigned char*)pbaTextures[2]+id,sizeof(unsigned char),cudaMemcpyDeviceToHost); //See which is the value at (sx,sy)
+
+    int iter = skelft2DFill(0,sx,sy,0,0,pbaTexSize,pbaTexSize,fill_value);                                //First, fill the background surrounding the image with some special value
+        
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+    
+
+    cudaBindTexture(0, pbaTexGray, pbaTextures[2]);                                                        //Used to read the thresholded DT
+
+    kernelFillHoles<<<grid, block>>>((unsigned char*)pbaTextures[2],pbaTexSize,background,foreground,fill_value);
+
+    cudaUnbindTexture(pbaTexGray);
+    
+    if (outputFill) cudaMemcpy(outputFill, (unsigned char*)pbaTextures[2], pbaTexSize * pbaTexSize * sizeof(unsigned char), cudaMemcpyDeviceToHost);
+    
+    return iter;
+}
+
+
+
+//--------  Advection-related code  ----------------------------------------------------------------------
+
+
+void advect2DSetSites(float* sites, int nsites)
+{
+    cudaMemcpy(pbaTexAdvectSites, sites, nsites * sizeof(float2), cudaMemcpyHostToDevice);                //Copy side-coords from CPU->GPU
+}
+
+//#define KERNEL(radius,dist) expf(-10*dist/radius)
+#define KERNEL(radius,dist) ((radius-dist)/radius)
+
+
+
+__global__ void kernelSplat(float* output, int size, int numpts, float d)
+//REMARK: This could be done tens of times faster using OpenGL (draw point-sprites with a 2D distance-kernel texture)
+{
+    int offs = blockIdx.x * blockDim.x + threadIdx.x;
+
+    if (offs < numpts)                                                        //careful not to index outside the site-vector..
+    {    
+      float2        p  = tex1Dfetch(pbaTexAdvect,offs);                        //splat the site whose coords are at location offs in pbaTexAdvect[]
+
+      float d2   = d*d;
+      short tx   = p.x, ty = p.y;                                            //REMARK: if I make these float, then something gets ccrewed up - why?
+      short jmin = ty-d;
+      short jmax = ty+d;
+
+      for(short j=jmin;j<=jmax;++j)                                            //bounding-box of the splat at 'skel'
+        {
+           float dy = (j-ty)*(j-ty);                                        //precalc this for speed
+           short imin = tx-d;
+           short imax = tx+d;
+           
+           int offs = __mul24(int(j),size);
+           
+           for(short i=imin;i<=imax;++i)                                    //REMARK: this could be made ~1.5x faster by computing the x-scanline intersection (min,max) with the d2-circle..
+             {
+                float r2 = dy+(i-tx)*(i-tx);
+                if (r2<=d2)                                                    //check we're really inside the splat
+                {
+                  int id = int(i)+offs; 
+                  float val = KERNEL(d2,r2);
+                  output[id] += val;                                        //Accumulate density (WARNING: hope this is thread-safe!!)
+                                                                            //If not, I need atomicAdd() on floats which requires CUDA 2.0    
+                }
+             }
+        }       
+    }
+}
+
+
+inline __device__ float2 computeGradient(const float2& p)                    //Compute -gradient of density at p
+{
+    float v_d = tex2D(pbaTexDensity,p.x,p.y-1);
+    float v_l = tex2D(pbaTexDensity,p.x-1,p.y);
+    float v_r = tex2D(pbaTexDensity,p.x+1,p.y);
+    float v_t = tex2D(pbaTexDensity,p.x,p.y+1);
+    
+    return make_float2((v_l-v_r)/2,(v_d-v_t)/2);
+}
+
+
+inline __device__ float2 computeGradientExplicit(float* density, const float2& p, int size)        //Compute gradient of density at p
+{
+    int x = p.x, y = p.y;
+    float v_d = density[TOID(x,y-1,size)];
+    float v_l = density[TOID(x-1,y,size)];
+    float v_r = density[TOID(x+1,y,size)];
+    float v_t = density[TOID(x,y+1,size)];
+    
+    return make_float2((v_r-v_l)/2,(v_t-v_d)/2);
+}
+
+
+__global__ void kernelGradient(float2* grad, int size)                            
+{
+    int tx = blockIdx.x * blockDim.x + threadIdx.x;
+    int ty = blockIdx.y * blockDim.y + threadIdx.y;
+
+    const float eps = 0.0001;
+
+    if (tx<size && ty<size)                                                    //Careful not to go outside the image..
+    {
+      float2 p = make_float2(tx,ty);
+      float2 g = computeGradient(p);
+      
+      float gn = sqrtf(g.x*g.x+g.y*g.y);                                    //robustly normalize the gradient
+      g.x /= gn+eps; g.y /= gn+eps;
+            
+      grad[TOID(tx,ty,size)] = g;    
+    }
+}
+
+
+void advect2DSplat(float* output, int num_pts, float radius)
+{
+    dim3 block = dim3(BLOCKX*BLOCKY);                                            //Prepare the splatting kernel: this operates on a vector of 2D sites
+    int numpts_b = (num_pts/block.x+1)*block.x;                                    //Find higher multiple of blocksize than # sites
+    dim3 grid  = dim3(numpts_b/block.x);
+    
+    cudaMemset(pbaTexDensityData,0,sizeof(float)*pbaTexSize*pbaTexSize);        //Zero the density texture
+    cudaBindTexture(0,pbaTexAdvect,pbaTexAdvectSites);                            //Bind sites to a texture
+
+    kernelSplat<<< grid, block >>>(pbaTexDensityData, pbaTexSize, num_pts, radius);
+                                                                                //Splat kernel    
+    cudaThreadSynchronize();    
+                                                                                
+    cudaUnbindTexture(pbaTexAdvect);                                            //Done with the sites
+                                                                                //Copy splat-map to CPU (if desired)
+    if (output) cudaMemcpy(output, pbaTexDensityData, pbaTexSize * pbaTexSize * sizeof(float), cudaMemcpyDeviceToHost);
+}
+
+
+
+void advect2DGetSites(float* output, int num_pts)                                //Get sites CUDA->CPU
+{
+    cudaMemcpy(output,pbaTexAdvectSites,num_pts*sizeof(float2),cudaMemcpyDeviceToHost);
+}
+
+
+void advect2DGradient()                                                            //Compute and save gradient of density map
+{
+    dim3 block = dim3(BLOCKX,BLOCKY);
+    dim3 grid  = dim3(pbaTexSize/block.x,pbaTexSize/block.y);
+    
+    cudaMemset(pbaTexGradientData,0,sizeof(float2)*pbaTexSize*pbaTexSize);        //Zero the gradient texture
+
+    pbaTexDensity.filterMode = cudaFilterModeLinear;                            //The floating-point density map
+    pbaTexDensity.normalized = false;
+    cudaChannelFormatDesc channelDesc1 = cudaCreateChannelDesc(32,0,0,0,cudaChannelFormatKindFloat);
+    cudaBindTexture2D(0,pbaTexDensity,pbaTexDensityData,channelDesc1,pbaTexSize,pbaTexSize,4*pbaTexSize);                            
+        
+    kernelGradient<<< grid, block >>>(pbaTexGradientData, pbaTexSize);            //Gradient kernel    
+                                                                                
+    cudaUnbindTexture(pbaTexDensity);
+}
+
+
+
+
+
-- 
cgit v1.2.3