API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Apr 6th, 2017
131
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 15.98 KB | None | 0 0
raw download clone embed print report
  1. /*
  2.  * Copyright 2011 Juan Gabriel Gomila Salas <joan-by@hotmail.com>
  3.  *
  4.  * All rights reserved.
  5.  *
  6.  *
  7.  * Patent warning:
  8.  *
  9.  * This file implement one algorithm possibly linked to the patent
  10.  *
  11.  * # N. Moroney, R.G. Beausoleil and I. Sobel, Local Color Correction,
  12.  * US Patent 6,822,762, (November 23, 2004)
  13.  *
  14.  * This file is made available for the exclusive aim of serving as
  15.  * scientific tool to verify the soundness and completeness of the
  16.  * algorithm description. Compilation, execution and redistribution
  17.  * of this file may violate patents rights in certain countries.
  18.  * The situation being different for every country and changing
  19.  * over time, it is your responsibility to determine which patent
  20.  * rights restrictions apply to you before you compile, use,
  21.  * modify, or redistribute this file. A patent lawyer is qualified
  22.  * to make this determination.
  23.  * If and only if they don't conflict with any patent terms, you
  24.  * can benefit from the following license terms attached to this
  25.  * file.
  26.  *
  27.  * License:
  28.  *
  29.  * This program is provided for scientific and educational only:
  30.  * you can use and/or modify it for these purposes, but you are
  31.  * not allowed to redistribute this work or derivative works in
  32.  * source or executable form. A license must be obtained from the
  33.  * patent right holders for any other use.
  34.  *
  35.  *
  36.  */
  37.  
  38. /**
  39.  * @file localcolorcorrection.c
  40.  * @brief Local Color Correction algorithm
  41.  *
  42.  * Parameters:
  43.  *
  44.  * input image
  45.  * output image
  46.  *
  47.  * Read/write operations (png format) make use
  48.  * of io_png.c and io_png.h, by Nicolas Limare
  49.  *
  50.  * Conversion operations between color spaces, make use
  51.  * colorspace library, by Pascal Getreuer
  52.  *
  53.  * @author Juan Gabriel Gomila Salas <joan-by@hotmail.com>
  54.  *
  55.  */
  56.  
  57. /**
  58.  * @mainpage Local Color Correction
  59.  *
  60.  * README.txt:
  61.  * @verbinclude README.txt
  62.  */
  63.  
  64. /**
  65.  * @file   localcolorcorrection.c
  66.  * @brief  Main executable file
  67.  *
  68.  *
  69.  * @author Juan Gabriel Gomila Salas <joan-by@hotmail.com>
  70.  */
  71.  
  72. #include <stdio.h>
  73. #include <stdlib.h>
  74. #include <string.h>
  75. #include <math.h>
  76. #include "structs.h"
  77. #include "io_png.h"
  78. #include <limits.h>
  79. #include <xmmintrin.h>  // SSE
  80. #include <pmmintrin.h>  // SSE
  81. #include <assert.h>
  82. #include <omp.h>
  83.  
  84.  
  85. #define MAX(X, Y) (((X) > (Y)) ? (X) : (Y))
  86. #define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
  87. /**
  88.  * @brief main function call
  89.  */
  90. void *contrast_enhance ( pixel *input, int w, int h )
  91. {
  92.     int   radius, n, nx, ny, t;
  93.     float *Mask, *Mask2, *Kernel;
  94.     float sum;
  95.  
  96.     pixel * output = (pixel *) malloc ( sizeof ( pixel ) * w * h );
  97.     radius = 40;
  98.  
  99.     //I     = (float *) malloc ( w * h * sizeof ( float ) );
  100.     Mask  = (float *) malloc ( w * h * sizeof ( float ) );
  101.     Mask2 = (float *) malloc ( w * h * sizeof ( float ) );
  102.     /*Create convolution kernel*/
  103.     Kernel = (float *) malloc ( ( 2 * radius + 1 ) * sizeof ( float ) );
  104. /*
  105.     for ( n = 0; n < w * h; n++ ) {
  106.         I[n] = ( ( input[n].R / 255.0 ) + ( input[n].G / 255.0 ) + ( input[n].B / 255.0 ) ) / 3.0f;
  107.     }
  108.  
  109.     for ( n = 0; n < w * h; n++ ) {
  110.         Mask[n]  = I[n];
  111.         Mask2[n] = Mask[n];
  112.     }*/
  113.  
  114.     /*#pragma omp parallel private(step,start,cpu,stop) shared(w,h,input,Mask,Mask2)
  115.     {
  116.       //step = (w*h) / omp_get_num_threads ();
  117.       //cpu = omp_get_thread_num ();
  118.       //stop = min ( ( cpu + 1 ) * step, (w*h));
  119.       //start = cpu * step;
  120.       #pragma omp for private(n)
  121.       for ( n = 0; n < w*h; n++) {
  122.         Mask[n] = ( ( input[n].R / 255.0 ) + ( input[n].G / 255.0 ) + ( input[n].B / 255.0 ) ) / 3.0f;
  123.         Mask2[n]  = Mask[n];
  124.     }
  125.   }*/
  126.   #pragma omp parallel for private(n) shared(w,h,input,Mask,Mask2)
  127.     for ( n = 0; n < w*h; n++) {
  128.       Mask[n] = ( ( input[n].R / 255.0 ) + ( input[n].G / 255.0 ) + ( input[n].B / 255.0 ) ) / 3.0f;
  129.       Mask2[n]  = Mask[n];
  130.   }
  131.  
  132.  
  133.   int step, cpu, stop, start;
  134. /*  #pragma omp parallel default(none) shared(input, Mask, Mask2, h, w) private(start, stop, cpu, step, n)
  135.   {
  136.     step = (w*h) / omp_get_num_threads();
  137.     cpu = omp_get_thread_num();
  138.     printf("Thread number %i\n", cpu);
  139.     stop = MIN ((cpu+1)*step, w*h);
  140.     start = cpu*step;
  141.       for ( n = start; n < stop; n++) {
  142.         Mask[n] = ( ( input[n].R / 255.0 ) + ( input[n].G / 255.0 ) + ( input[n].B / 255.0 ) ) / 3.0f;
  143.         Mask2[n]  = Mask[n];
  144.     }
  145.   }*/
  146.  
  147.     if ( radius > 0 ) {
  148.         #pragma omp parallel for default(none) private(n) shared(radius, Kernel)
  149.         for ( n = 0; n < 2 * radius + 1; n++ ) {
  150.             Kernel[n] =
  151.                 (float) exp ( (double) ( -2.0f * ( -radius + n ) * ( -radius + n ) / ( radius * radius ) ) );
  152.         }
  153.  
  154.         sum = 0.0f;
  155.         //#pragma omp parallel for default(none) reduction(+:sum) shared(radius, Kernel)
  156.         for ( n = 0; n < 2 * radius + 1; n++ ) {
  157.                    sum = sum + Kernel[n];
  158.         }
  159.         #pragma omp parallel for default(none) private(n) shared(radius, Kernel, sum)
  160.         for ( n = 0; n < 2 * radius + 1; n++ ) {
  161.             Kernel[n] = (float) ( Kernel[n] / sum );
  162.         }
  163.  
  164.         /*If the radius of the mask is less than the mid-size of the image,
  165.            we perform the desired algorithm*/
  166.         if ( radius < w / 2 && radius < h / 2 ) {
  167.           /*Convolve in the x direction*/
  168.         /*  #pragma omp parallel for default(none) shared(radius,w,h,Mask,Kernel,Mask2) private(ny,t, sum)
  169.             for ( nx = 0; nx < w; nx++ ) {
  170.                 for ( ny = 0; ny < h; ny++ ) {
  171.                     sum = 0.0f;
  172.                     for ( t = -radius; t <= radius; t++ ) {
  173.                         if ( nx + t >= w ) {
  174.                             sum = sum + Mask[2 * w - 2 - nx - t + ny * w] * Kernel[t + radius];
  175.                         }
  176.                         else{
  177.                             sum = sum + Mask[abs ( nx + t ) + ny * w] * Kernel[t + radius];
  178.                         }
  179.                     }
  180.                     Mask2[nx + ny * w] = (float) sum;
  181.                 }
  182.             }*/
  183.  
  184. /*
  185.     //WITHOUT ABS()
  186.  #define MAX(X, Y) (((X) > (Y)) ? (X) : (Y))
  187.  #define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
  188.      int limit, limit2;
  189.      for ( nx = 0; nx < w; nx++ ) {
  190.          for ( ny = 0; ny < h; ny++ ) {
  191.              sum = 0.0f;
  192.              for ( t = w-nx; t <= radius; t++ ) {
  193.                 sum = sum + Mask[2 * w - 2 - nx - t + ny * w] * Kernel[t + radius];
  194.              }
  195.              limit2 = MAX(-nx, -radius-1);
  196.              for ( t = -radius; t <= limit2; t++ ) {
  197.                  sum = sum + Mask[ -(nx + t ) + ny * w] * Kernel[t + radius];
  198.              }
  199.              limit = MIN(w-nx-1, radius);
  200.              for ( t = limit2+1; t <= limit; t++ ) {
  201.                  sum = sum + Mask[( nx + t ) + ny * w] * Kernel[t + radius];
  202.              }
  203.              Mask2[nx + ny * w] = (float) sum;
  204.          }
  205.      }
  206. */
  207.      //WITH OpenMP
  208.     /* int limit, limit2;
  209.      #pragma omp parallel for private(t, ny, limit, limit2, sum) shared(w,h,Mask,Mask2,Kernel, radius) //collapse(2)
  210.      for ( nx = 0; nx < w; nx++ ) {
  211.        for ( ny = 0; ny < h; ny++ ) {
  212.          sum = 0.0f;
  213.          for ( t = w-nx; t <= radius; t++ ) {
  214.            sum = sum + Mask[2 * w - 2 - nx - t + ny * w] * Kernel[t + radius];
  215.          }
  216.          limit2 = MAX(-nx, -radius-1);
  217.          for ( t = -radius; t <= limit2; t++ ) {
  218.            sum = sum + Mask[ -(nx + t ) + ny * w] * Kernel[t + radius];
  219.          }
  220.          limit = MIN(w-nx-1, radius);
  221.          for ( t = limit2+1; t <= limit; t++ ) {
  222.            sum = sum + Mask[( nx + t ) + ny * w] * Kernel[t + radius];
  223.          }
  224.          Mask2[nx + ny * w] = (float) sum;
  225.        }
  226.      }*/
  227.  
  228.      int limit, limit2;
  229.      //#pragma omp parallel private(t, ny, limit, limit2, sum) shared(w,h,Mask,Mask2,Kernel, radius) //collapse(2)
  230.      #pragma omp parallel private(t, ny, nx, limit, limit2, sum, start,step, cpu, stop) shared(w,h,Mask,Mask2,Kernel,radius) //collapse(2)
  231.      {
  232.        step = (w) / omp_get_num_threads();
  233.        cpu = omp_get_thread_num();
  234.        stop = MIN ((cpu+1)*step, w);
  235.        start = cpu*step;
  236.        for ( nx = start; nx < stop; nx++ ) {
  237.          for ( ny = 0; ny < h; ny++ ) {
  238.            sum = 0.0f;
  239.            for ( t = w-nx; t <= radius; t++ ) {
  240.              sum = sum + Mask[2 * w - 2 - nx - t + ny * w] * Kernel[t + radius];
  241.            }
  242.            limit2 = MAX(-nx, -radius-1);
  243.            for ( t = -radius; t <= limit2; t++ ) {
  244.              sum = sum + Mask[ -(nx + t ) + ny * w] * Kernel[t + radius];
  245.            }
  246.            limit = MIN(w-nx-1, radius);
  247.            for ( t = limit2+1; t <= limit; t++ ) {
  248.              sum = sum + Mask[( nx + t ) + ny * w] * Kernel[t + radius];
  249.            }
  250.            Mask2[nx + ny * w] = (float) sum;
  251.          }
  252.        }
  253.      }
  254.      /*int limit, limit2;
  255.      #pragma omp parallel for private(t, ny, limit, limit2) shared(w,h,Mask,Mask2,Kernel, radius,start, stop, cpu, step) //collapse(2)
  256.      for ( nx = 0; nx < w; nx++ ) {
  257.      #pragma omp parallel default(none) shared(Kernel, Mask, Mask2, h, w, nx, radius) private(start, stop, cpu, step, ny, sum,t, limit, limit2)
  258.      {
  259.        step = (h) / omp_get_num_threads();
  260.        cpu = omp_get_thread_num();
  261.        printf("Thread number %i\n", cpu);
  262.        stop = MIN ((cpu+1)*step, h);
  263.        start = cpu*step;
  264.        for ( ny = start; ny < stop; ny++ ) {
  265.          sum = 0.0f;
  266.          for ( t = w-nx; t <= radius; t++ ) {
  267.            sum = sum + Mask[2 * w - 2 - nx - t + ny * w] * Kernel[t + radius];
  268.          }
  269.          limit2 = MAX(-nx, -radius-1);
  270.          for ( t = -radius; t <= limit2; t++ ) {
  271.            sum = sum + Mask[ -(nx + t ) + ny * w] * Kernel[t + radius];
  272.          }
  273.          limit = MIN(w-nx-1, radius);
  274.          for ( t = limit2+1; t <= limit; t++ ) {
  275.            sum = sum + Mask[( nx + t ) + ny * w] * Kernel[t + radius];
  276.          }
  277.          Mask2[nx + ny * w] = (float) sum;
  278.        }
  279.       }
  280.     }*/
  281.  
  282.  
  283. /*#define MAX(X, Y) (((X) > (Y)) ? (X) : (Y))
  284. #define MIN(X, Y) (((X) < (Y)) ? (X) : (Y))
  285. int limit, limit2;
  286. #pragma omp parallel for default(none) shared(w) private(nx)
  287. for ( nx = 0; nx < w; nx++ ) {
  288.   #pragma omp parallel for default(none) shared(h, w, radius, Mask2, nx) private(ny, sum, limit, limit2)
  289.   for ( ny = 0; ny < h; ny++ ) {
  290.     sum = 0.0f;
  291.     #pragma omp parallel for reduction(+:sum) default(none) shared(radius, nx, ny, Mask, Kernel, w) private(t)
  292.     for ( t = w-nx; t <= radius; t++ ) {
  293.       sum = sum + Mask[2 * w - 2 - nx - t + ny * w] * Kernel[t + radius];
  294.     }
  295.     limit2 = MAX(-nx, -radius-1);
  296.     #pragma omp parallel for reduction(+:sum) default(none) shared(radius, nx, ny, Mask, Kernel, limit2, w) private(t)
  297.     for ( t = -radius; t <= limit2; t++ ) {
  298.       sum = sum + Mask[ -(nx + t ) + ny * w] * Kernel[t + radius];
  299.     }
  300.     limit = MIN(w-nx-1, radius);
  301.     #pragma omp parallel for reduction(+:sum) default(none) shared(radius, nx, ny, Mask, Kernel, limit2, limit, w) private(t)
  302.     for ( t = limit2+1; t <= limit; t++ ) {
  303.       sum = sum + Mask[( nx + t ) + ny * w] * Kernel[t + radius];
  304.     }
  305.     Mask2[nx + ny * w] = (float) sum;
  306.   }
  307. }/*
  308.             /*Convolve in the y direction*/
  309.             /*#pragma omp parallel for default(none) shared(radius,w,h,Mask,Kernel,Mask2) private(nx,t, sum)
  310.             for ( ny = 0; ny < h; ny++ ) {
  311.                 for ( nx = 0; nx < w; nx++ ) {
  312.                     sum = 0.0f;
  313.                     for ( t = -radius; t <= radius; t++ ) {
  314.                         if ( ny + t >= h ) {
  315.                             sum = sum + Mask2[nx + ( 2 * h - 2 - ny - t ) * w] *
  316.                                   Kernel[t + radius];
  317.                         }
  318.                         else{
  319.                             sum = sum + Mask2[nx + abs ( ny + t ) * w] *
  320.                                   Kernel[t + radius];
  321.                         }
  322.                     }
  323.                     Mask[nx + ny * w] = (float) sum;
  324.                 }
  325.             }*/
  326.             /*#pragma omp parallel for default(none) private(t, nx, limit, limit2, sum) shared(w,h,Mask,Mask2,Kernel, radius) // collapse(2)
  327.             for ( ny = 0; ny < h; ny++ ) {
  328.                 for ( nx = 0; nx < w; nx++ ) {
  329.                     sum = 0.0f;
  330.                     for ( t = h-ny; t <= radius; t++ ) {
  331.                       sum = sum + Mask2[nx + ( 2 * h - 2 - ny - t ) * w] * Kernel[t + radius];
  332.                     }
  333.                     limit2 = MAX(-ny, -radius-1);
  334.                     for ( t = -radius; t <= limit2; t++ ) {
  335.                       sum = sum + Mask2[nx - ( ny + t ) * w] * Kernel[t + radius];
  336.                     }
  337.                     limit = MIN(h-ny-1, radius);
  338.                     for ( t = limit2+1; t <= limit; t++ ){
  339.                       sum = sum + Mask2[nx + ( ny + t ) * w] * Kernel[t + radius];
  340.                     }
  341.                     Mask[nx + ny * w] = (float) sum;
  342.                 }
  343.             }*/
  344.  
  345.             #pragma omp parallel private(t, ny, nx, limit, limit2, sum, start,step, cpu, stop) shared(w,h,Mask,Mask2,Kernel,radius) //collapse(2)
  346.             {
  347.               step = (h) / omp_get_num_threads();
  348.               cpu = omp_get_thread_num();
  349.               stop = MIN ((cpu+1)*step, h);
  350.               start = cpu*step;
  351.               for ( ny = start; ny < stop; ny++ ) {
  352.                   for ( nx = 0; nx < w; nx++ ) {
  353.                       sum = 0.0f;
  354.                       for ( t = h-ny; t <= radius; t++ ) {
  355.                         sum = sum + Mask2[nx + ( 2 * h - 2 - ny - t ) * w] * Kernel[t + radius];
  356.                       }
  357.                       limit2 = MAX(-ny, -radius-1);
  358.                       for ( t = -radius; t <= limit2; t++ ) {
  359.                         sum = sum + Mask2[nx - ( ny + t ) * w] * Kernel[t + radius];
  360.                       }
  361.                       limit = MIN(h-ny-1, radius);
  362.                       for ( t = limit2+1; t <= limit; t++ ){
  363.                         sum = sum + Mask2[nx + ( ny + t ) * w] * Kernel[t + radius];
  364.                       }
  365.                       Mask[nx + ny * w] = (float) sum;
  366.                   }
  367.               }
  368.             }
  369.  
  370.         }
  371.         else {
  372.             /*Otherwise, perform a simple gamma correction*/
  373.             //#pragma omp parallel for default(none) private(nx) shared(Mask, h, w) reduction(+:sum) collapse(2)
  374.             for ( ny = 0; ny < h; ny++ ) {
  375.                 for ( nx = 0; nx < w; nx++ ) {
  376.                     sum = sum + Mask[nx + ny * w];
  377.                 }
  378.             }
  379.             sum = sum / ( w * h );
  380.             #pragma omp parallel for default(none) private(nx) shared(Mask, h, w, sum)  collapse(2)
  381.             for ( ny = 0; ny < h; ny++ ) {
  382.                 for ( nx = 0; nx < w; nx++ ) {
  383.                     Mask[nx + ny * w] = sum;
  384.                 }
  385.             }
  386.         }
  387.     }
  388.  
  389.   /*  #pragma omp parallel for default(none) private(n) shared(Mask, output, input, h, w)
  390.     for ( n = 0; n < w * h; n++ ) {
  391.         Mask[n] = (float) ( 2 * Mask[n] - 1 );
  392.         Mask[n] = (float) pow ( 2.0, (double) Mask[n] );
  393.  
  394.         output[n].R = min ( 255, 255 * (float) pow ( (double) input[n].R / 255.0, (double) Mask[n] ) );
  395.         output[n].G = min ( 255, 255 * (float) pow ( (double) ( input[n].G / 255.0 ), (double) Mask[n] ) );
  396.         output[n].B = min ( 255, 255 * (float) pow ( (double) ( input[n].B / 255.0 ), (double) Mask[n] ) );
  397.     }*/
  398.     #pragma omp parallel default(none) private(n,start,step, cpu, stop) shared(Mask, output, input, h, w)
  399.     {
  400.       step = (h*w) / omp_get_num_threads();
  401.       cpu = omp_get_thread_num();
  402.       stop = MIN ((cpu+1)*step, h*w);
  403.       start = cpu*step;
  404.       for ( n = start; n < stop; n++ ) {
  405.           Mask[n] = (float) ( 2 * Mask[n] - 1 );
  406.           Mask[n] = (float) pow ( 2.0, (double) Mask[n] );
  407.  
  408.           output[n].R = min ( 255, 255 * (float) pow ( (double) input[n].R / 255.0, (double) Mask[n] ) );
  409.           output[n].G = min ( 255, 255 * (float) pow ( (double) ( input[n].G / 255.0 ), (double) Mask[n] ) );
  410.           output[n].B = min ( 255, 255 * (float) pow ( (double) ( input[n].B / 255.0 ), (double) Mask[n] ) );
  411.       }
  412.     }
  413.     //free ( I );
  414.     free ( Mask );
  415.     free ( Mask2 );
  416.     free ( Kernel );
  417.  
  418.     return output;
  419. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!