source: proiecte/pmake3d/segment/segment.cu @ 77

// example1.cpp : Defines the entry point for the console application.
//

//#include "stdafx.h"

#include <vector>
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <cuda.h>
#include <cudpp.h>
#include <cutil.h>
#include <math_functions.h>
#include <image.h>
#include <misc.h>
#include <pnmfile.h>
#include <cutil_math.h>
#include <sys/time.h>
#include <unistd.h>


#include "segment-image.h"


#define BLOCK_DIM_X 12
#define BLOCK_DIM_Y 8

#define max FLT_MAX



__global__ void convolve_even_gpu(imagef *src, imagef *dst,  float *mask, int dim_mask){
        int w = src->w;
        int h = src->h;
        int x = blockIdx.x * blockDim.x + threadIdx.x;
        int y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx = x + y*w;
        int i;
        float sum;      
        int mmax,mmin;

        if (x>=w || y>=h) 
                return;
        
        sum = mask[0] * src->data[idx];

        for (i = 1; i < dim_mask; i++) {
                mmax = (x-i > 0 ? x-i : 0);
                //mmax = max(x-i,0);// > 0 ? x-i : 0)
                mmin = (x+i > w-1 ? x+i : w-1);
                sum += mask[i] * (src->data[mmax + y * w]) + src->data[mmin + y * w];
    }
    dst->data[y + h * x] = sum;

}


__global__ void create_rgb_smooth(imagef *r, imagef *g, imagef *b, image *input)
{
        int w = input->w;
        int h = input->h;
        int x = blockIdx.x * blockDim.x + threadIdx.x;
        int y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx = x + y*w;
        
        if (x>=w || y>=h) 
                return; 

        r->data[idx] = input->data[idx].r;
        g->data[idx] = input->data[idx].g;
        b->data[idx] = input->data[idx].b;

}

#define diffgpu(r, g, b, x1, y1, x2, y2)                                        \
  sqrtf(powf(r->data[x1+w*y1]-r->data[x2+w*y2],2) +                     \
              powf(g->data[x1+w*y1]-g->data[x2+w*y2],2) +               \
              powf(b->data[x1+w*y1]-b->data[x2+w*y2],2));               \


// Kernel that executes on the CUDA device
__global__ void compute_edges(image *input, edge *edges, edge **edges_location,float *weights,imagef *smooth_r,
                                        imagef *smooth_g,imagef *smooth_b )
{
        //int idx = blockIdx.x * blockDim.x + threadIdx.x;
        int w = input->w;
        int h = input->h;
        int x = blockIdx.x * blockDim.x + threadIdx.x;
        int y = blockIdx.y * blockDim.y + threadIdx.y;
        int idx = x + y*w;

//      indx[idx] = idx;

/*      if (blockIdx.x == 0 && blockIdx.y == 0 && threadIdx.x == 0 && threadIdx.y == 0)
                input->w = input->w + 1;

        edges[0].a = 23;
        edges[0].b = 24;
*/

        if (x>=w || y>=h) 
                return; 

        if (x < w-1) 
        {
                edges[idx*4].a = y * w + x;
                edges[idx*4].b = y * w + (x+1);         
                weights[idx*4] = diffgpu(smooth_r, smooth_g, smooth_b, x, y, x+1, y);
        }
        else
        {       
                edges[idx*4].a = -1;
                edges[idx*4].b = -1;
                weights[idx*4] = max;
        }
                

    if (y < h-1) 
        {
                edges[idx*4+1].a = y * w + x;
                edges[idx*4+1].b = (y+1) * w + x;
                weights[idx*4+1] = diffgpu(smooth_r, smooth_g, smooth_b, x, y, x, y+1);
    }
        else
        {
                edges[idx*4+1].a = -1;
                edges[idx*4+1].b = -1;
                weights[idx*4+1] = max;
        }

    if ((x < w-1) && (y < h-1)) 
        {
                edges[idx*4+2].a = y * w + x;
                edges[idx*4+2].b = (y+1) * w + (x+1);
                weights[idx*4+2] = diffgpu(smooth_r, smooth_g, smooth_b, x, y, x+1, y+1);
        }
        else
        {
                edges[idx*4+2].a = -1;
                edges[idx*4+2].b = -1;
                weights[idx*4+2] = max;
        
        }

    if ((x < w-1) && (y > 0)) 
        {
                edges[idx*4+3].a = y * w + x;
                edges[idx*4+3].b = (y-1) * w + (x+1);
                weights[idx*4+3] = diffgpu(smooth_r, smooth_g, smooth_b, x, y, x+1, y-1);
    }
        else
        {
                edges[idx*4+3].a = -1;
                edges[idx*4+3].b = -1;
                weights[idx*4+3] = max;
        }

        edges_location[idx*4] = &(edges[idx*4]);    
        edges_location[idx*4+1] = &(edges[idx*4+1]);    
        edges_location[idx*4+2] = &(edges[idx*4+2]);    
        edges_location[idx*4+3] = &(edges[idx*4+3]);    


        /*  if (idx == 0)
                        input->h = input->h + 1;
        */              

}

// main routine that executes on the host
int main(int argc, char **argv)
{
        if (argc != 6) 
        {
                fprintf(stderr, "usage: %s sigma k min input(ppm) output(ppm)\n", argv[0]);
                return 1;
        }
 
        
        dim3 n_blocks, block_size; 
        //int n_blocks, block_size; 
        float sigma = atof(argv[1]);
        float k = atof(argv[2]);
        int min_size = atoi(argv[3]);
        rgb *array;
        edge *edges;
        edge *ed_temp;

        float* weights; 
        edge** edges_location;

        
        struct timeval start, end;
        long mtime, seconds, useconds;    

    gettimeofday(&start, NULL);

        printf("loading input image.\n");
        image *input = loadPPM(argv[4]);

        /* Alloc image on device*/
        image *input_cuda;

        image in;

        cudaMalloc((void **) &array, input->h * input->w * sizeof(rgb));
        cudaMemcpy((void *) array, (void *) input->data, input->h * input->w * sizeof(rgb), cudaMemcpyHostToDevice);    
        

        in.h = input->h;
        in.w = input->w;
        in.data = (rgb*)array;
  
        cudaMalloc((void **) &input_cuda, sizeof(image));       
        cudaMemcpy((void *) input_cuda, (void *) &in, sizeof(image), cudaMemcpyHostToDevice);


        /* Alocare matrci r,g,b de image_float*/
        imagef *r; 
        imagef *g; 
        imagef *b;
//      imagef *tmp;
        imagef tmp;

        imagef *tmp_img;
        imagef tmp_img_cpu;
        
        float *arrayftmp;
        
        
        float *arrayf1, *arrayf2, *arrayf3;
        std::vector<float> vmask = make_fgauss(sigma);
        float *mask = (float *)malloc(vmask.size()*sizeof(float)) + 3 * sizeof(imagef);
        for (int i=0; i<vmask.size(); i++)
                mask[i]=vmask[i];
        float *maski;
//      tmp = (imagef *)(((char *)mask) + vmask.size()*sizeof(float));
        
                        

    unsigned int rgb_size = 4 * input->h * input->w * sizeof(float) + 4 * sizeof(imagef) + sizeof(float) * vmask.size();
        void *bigMem;   
        if (cudaMalloc(&bigMem, rgb_size) != cudaSuccess)
                printf("ERROR big_malloc 1\n");
        char *initial = (char *)bigMem;
        maski = (float *)initial;
        initial += sizeof(float) * vmask.size();        
        r = (imagef *)initial;
        initial+=sizeof(imagef);        
        g = (imagef *)initial;
        initial+=sizeof(imagef);
        b = (imagef *)initial;
        initial+=sizeof(imagef);
        tmp_img = (imagef *)initial;
        initial+=sizeof(imagef);
        arrayf1 = (float *)initial;
        initial += input->h * input->w * sizeof(float);
        arrayf2 = (float *)initial;
        initial += input->h * input->w * sizeof(float);
        arrayf3 = (float *)initial;
        initial += input->h * input->w * sizeof(float);
        arrayftmp = (float *)initial;
/*      
        for (int i=0; i<3; i++){
                tmp[i].w=input->w;
                tmp[i].h=input->h;      
        }
        tmp[0].data = arrayf1;
        tmp[1].data = arrayf2;
        tmp[2].data = arrayf3;
        CUDA_SAFE_CALL(cudaMemcpy((void *)maski, (void *)mask, 3*sizeof(imagef) + vmask.size()*sizeof(float), cudaMemcpyHostToDevice));
*/      
//      exit(0);
/*
    if (cudaMalloc((void **) &r, sizeof(imagef)) != cudaSuccess)
                printf("ERROR11\n");
        if (cudaMalloc((void **) &g, sizeof(imagef)) != cudaSuccess)
                printf("ERROR12\n");

        if (cudaMalloc((void **) &b, sizeof(imagef)) != cudaSuccess)
                printf("ERROR13\n");

        cudaMalloc((void **) &arrayf1, input->h * input->w * sizeof(float));
*/      

        tmp.w = input->w;
        tmp.h = input->h;
        tmp.data = arrayf1;

        CUDA_SAFE_CALL(cudaMemcpy((void *) r, (void *) &tmp, sizeof(imagef), cudaMemcpyHostToDevice));  
        
//      cudaMalloc((void **) &arrayf2, input->h * input->w * sizeof(float));
        tmp.data = arrayf2;
        if (cudaMemcpy((void *) g, (void *) &tmp, sizeof(imagef), cudaMemcpyHostToDevice)!=cudaSuccess)
                printf("Error memcpy g\n");     
        
//      cudaMalloc((void **) &arrayf3, input->h * input->w * sizeof(float));
        tmp.data = arrayf3;
        if (cudaMemcpy((void *) b, (void *) &tmp, sizeof(imagef), cudaMemcpyHostToDevice)!=cudaSuccess)
                printf("Error memcpy b\n");     
        
        if (cudaMemcpy((void *) maski, (void *) mask, vmask.size() * sizeof(float), cudaMemcpyHostToDevice) != cudaSuccess)
                printf("Error memcpy maski\n");

        
        //smooth alloc

        imagef *smooth_r; 
        imagef *smooth_g; 
        imagef *smooth_b;
        imagef smooth_tmp;
        
        
        float *arrayfs1, *arrayfs2, *arrayfs3;

        
        unsigned int smooth_size = 3 * input->h * input->w * sizeof(float) + 3 * sizeof(imagef);
        void *bigMemSmooth;     
        CUDA_SAFE_CALL(cudaMalloc(&bigMemSmooth, smooth_size));

        initial = (char *)bigMemSmooth;
        smooth_r = (imagef *)initial;
        initial+=sizeof(imagef);        
        smooth_g = (imagef *)initial;
        initial+=sizeof(imagef);
        smooth_b = (imagef *)initial;
        initial+=sizeof(imagef);
        arrayfs1 = (float *)initial;
        initial += input->h * input->w * sizeof(float);
        arrayfs2 = (float *)initial;
        initial += input->h * input->w * sizeof(float);
        arrayfs3 = (float *)initial;
        
/*
    if (cudaMalloc((void **) &smooth_r, sizeof(imagef)) != cudaSuccess)
                printf("ERROR14\n");
        if (cudaMalloc((void **) &smooth_g, sizeof(imagef)) != cudaSuccess)
                printf("ERROR15\n");

        if (cudaMalloc((void **) &smooth_b, sizeof(imagef)) != cudaSuccess)
                printf("ERROR16\n");

        if (cudaMalloc((void **) &arrayfs1, input->h * input->w * sizeof(float))!=cudaSuccess)
                printf("Error malloc arrayfs1\n");
*/
        smooth_tmp.w = input->w;
        smooth_tmp.h = input->h;
        smooth_tmp.data = arrayfs1;
        if (cudaMemcpy((void *) smooth_r, (void *) &smooth_tmp, sizeof(imagef), cudaMemcpyHostToDevice)!=cudaSuccess)
                printf("Error memcpy smoothr\n");       
        
//      if (cudaMalloc((void **) &arrayfs2, input->h * input->w * sizeof(float))!=cudaSuccess)
//              printf("Error malloc arrayfs2\n");
        smooth_tmp.data = arrayfs2;
        if (cudaMemcpy((void *) smooth_g, (void *) &smooth_tmp, sizeof(imagef), cudaMemcpyHostToDevice)!=cudaSuccess)
                printf("Error memcpy smoothg\n");       
        
//      if (cudaMalloc((void **) &arrayfs3, input->h * input->w * sizeof(float))!=cudaSuccess)
//              printf("Error malloc arrayfs3\n");
        smooth_tmp.data = arrayfs3;
        if (cudaMemcpy((void *) smooth_b, (void *) &smooth_tmp, sizeof(imagef), cudaMemcpyHostToDevice)!=cudaSuccess)
                printf("Error memcpy smoothb\n");       
        
        //alloc tmp matrix for smooth
//      if (cudaMalloc((void **) &tmp_img, sizeof(imagef)) != cudaSuccess)
//              printf("ERROR17\n");
//      cudaMalloc((void **) &arrayftmp, input->w * input->h * sizeof(float));
        tmp_img_cpu.w = input->h;
        tmp_img_cpu.h = input->w;
        tmp_img_cpu.data = arrayftmp;
        if (cudaMemcpy((void *) tmp_img, (void *) &tmp_img_cpu, sizeof(imagef), cudaMemcpyHostToDevice)!=cudaSuccess)
                printf("Error memcpy tmp_img\n");       
        
//      if (cudaMalloc((void **) &maski, sizeof(float)*vmask.size()) != cudaSuccess)
//              printf("ERROR18\n");
        
        
        create_rgb_smooth <<< n_blocks, block_size >>> (r, g, b, input_cuda);
        
        convolve_even_gpu <<< n_blocks, block_size >>> (r, tmp_img, maski, vmask.size());
        convolve_even_gpu <<< n_blocks, block_size >>> (tmp_img, smooth_r, maski, vmask.size());
        convolve_even_gpu <<< n_blocks, block_size >>> (g, tmp_img, maski, vmask.size());
        convolve_even_gpu <<< n_blocks, block_size >>> (tmp_img, smooth_g, maski, vmask.size());
        convolve_even_gpu <<< n_blocks, block_size >>> (b, tmp_img, maski, vmask.size());
        convolve_even_gpu <<< n_blocks, block_size >>> (tmp_img, smooth_b, maski, vmask.size());
        

        //free(mask);
        /*      
        cudaFree(maski);
        cudaFree(tmp_img);
        cudaFree(r);
        cudaFree(g);
        cudaFree(b);
        cudaFree(arrayf1);
        cudaFree(arrayf2);
        cudaFree(arrayf3);
        cudaFree(arrayftmp);
        */
        CUDA_SAFE_CALL(cudaFree(input_cuda));   
        CUDA_SAFE_CALL(cudaFree(array));        
        
        CUDA_SAFE_CALL(cudaFree(bigMem));       
        
        
        /* Alloc edges' vector on device */
        //CUDA_SAFE_CALL(cudaMalloc((void **) &edges, (4 * input->h * input->w )* sizeof(edge)));       
        //CUDA_SAFE_CALL(cudaMalloc((void **) &weights, (4 * input->h * input->w )* sizeof(float)));
        unsigned int size_edges = 4 * input->h * input->w * sizeof(edge) + 4 * input->h * input->w * sizeof(edge*) + 4 * input->h * input->w * sizeof(float);
        CUDA_SAFE_CALL(cudaMalloc((void **) &bigMem, size_edges));
        
        //exit(0);

        /*
        if (cudaMalloc((void **) &edges, 4 * input->h * input->w * sizeof(edge)) != cudaSuccess)
                printf("ERROR19\n");    
        if (cudaMalloc((void **) &edges_location, (4 * input->h * input->w )* sizeof(edge*)) != cudaSuccess)    
                printf("ERROR1A\n");    
        
        if (cudaMalloc((void **) &weights, (4 * input->h * input->w )* sizeof(float)) != cudaSuccess)   
                printf("ERROR1B\n");    
        */
        initial = (char *)bigMem;
        edges = (edge*)initial;
        initial += 4 * input->h * input->w * sizeof(edge);
        edges_location = (edge **)initial;
        initial += 4 * input->h * input->w * sizeof(edge*);
        weights = (float *)initial;

    /* Define the grid */
        n_blocks.x = input->w / BLOCK_DIM_X + (input->w % BLOCK_DIM_X == 0?0:1); 
        n_blocks.y = input->h / BLOCK_DIM_Y + (input->h % BLOCK_DIM_Y == 0?0:1); 
        n_blocks.z = 1;

        printf("n_blocks.x = %d\n", n_blocks.x);
        printf("n_blocks.y = %d\n", n_blocks.y);


        block_size.x = BLOCK_DIM_X;
        block_size.y = BLOCK_DIM_Y;
        block_size.z = 1;
        
        /* Compute edges */
        compute_edges <<< n_blocks, block_size >>> (input_cuda, edges, edges_location, weights,smooth_r,smooth_g,smooth_b);
/*
        cudaFree(smooth_r);
        cudaFree(smooth_g);
        cudaFree(smooth_b);
        cudaFree(arrayfs1);
        cudaFree(arrayfs2);
        cudaFree(arrayfs3);
*/
        CUDA_SAFE_CALL(cudaFree(bigMemSmooth));
        
        ed_temp = (edge *)malloc((4 * input->h * input->w )* sizeof(edge));
        if (cudaMemcpy((void *) ed_temp, (void *) edges, 4 * input->h * input->w * sizeof(edge), cudaMemcpyDeviceToHost) != cudaSuccess)
                printf("ERROR2\n");

/*      if (cudaMemcpy((void *) indx_temp, (void *) indx, input->h * input->w * sizeof(int), cudaMemcpyDeviceToHost))
                printf("ERROR3\n");     
*/

        printf("a0 = %d\n", ed_temp[0].a);
        printf("b0 = %d\n", ed_temp[0].b);
        
        //sort edges
        /*
        CUDPPConfiguration config;
    config.op = CUDPP_MIN;
    config.datatype = CUDPP_FLOAT;
    config.algorithm = CUDPP_SORT_RADIX;
    config.options = CUDPP_OPTION_KEY_VALUE_PAIRS;
    
    CUDPPHandle scanplan = 0;
    CUDPPResult result = cudppPlan(&scanplan, config, 4 * input->h * input->w, 1, 0);  
        
        result = cudppSort (  scanplan,
                weights,
                edges_location,
                32,
                4 * input->h * input->w
        );
        
        result = cudppDestroyPlan(scanplan);
    if (CUDPP_SUCCESS != result)
    {
        printf("Error destroying CUDPPPlan\n");
        exit(-1);
    }


    if (CUDPP_SUCCESS != result)
    {
        printf("Error creating CUDPPPlan\n");
        exit(-1);
    }
        */
        gettimeofday(&end, NULL);

    seconds  = end.tv_sec  - start.tv_sec;
    useconds = end.tv_usec - start.tv_usec;

    mtime = ((seconds) * 1000 + useconds/1000.0) + 0.5;

    printf("Elapsed time: %ld milliseconds\n", mtime);
        
        // Cleanup
        //exit(0);
        //cudaFree(input_cuda); 
        CUDA_SAFE_CALL(cudaFree(bigMem));
        //cudaFree(array);
        //cudaFree(edges);
        //cudaFree(edges_location);
        //cudaFree(weights);    


        return 0;

}