src/display/nr-filter-turbulence.cpp

   1 /*
   2  * feTurbulence filter primitive renderer
   3  *
   4  * Authors:
   5  *   World Wide Web Consortium <http://www.w3.org/>
   6  *   Felipe Corrêa da Silva Sanches <juca@members.fsf.org>
   7  *
   8  * This file has a considerable amount of code adapted from
   9  *  the W3C SVG filter specs, available at:
  10  *  http://www.w3.org/TR/SVG11/filters.html#feTurbulence
  11  *
  12  * W3C original code is licensed under the terms of
  13  *  the (GPL compatible) W3C® SOFTWARE NOTICE AND LICENSE:
  14  *  http://www.w3.org/Consortium/Legal/2002/copyright-software-20021231
  15  *
  16  * Copyright (C) 2007 authors
  17  * Released under GNU GPL version 2 (or later), read the file 'COPYING' for more information
  18  */
  19
  20 #include "display/nr-arena-item.h"
  21 #include "display/nr-filter.h"
  22 #include "display/nr-filter-turbulence.h"
  23 #include "display/nr-filter-units.h"
  24 #include "display/nr-filter-utils.h"
  25 #include "libnr/nr-rect-l.h"
  26 #include "libnr/nr-blit.h"
  27 #include <math.h>
  28
  29 namespace Inkscape {
  30 namespace Filters{
  31
  32 FilterTurbulence::FilterTurbulence()
  33 : XbaseFrequency(0),
  34   YbaseFrequency(0),
  35   numOctaves(1),
  36   seed(0),
  37   updated(false),
  38   updated_area(NR::IPoint(), NR::IPoint()),
  39   pix(NULL),
  40   fTileWidth(10), //guessed
  41   fTileHeight(10), //guessed
  42   fTileX(1), //guessed
  43   fTileY(1) //guessed
  44 {
  45 }
  46
  47 FilterPrimitive * FilterTurbulence::create() {
  48     return new FilterTurbulence();
  49 }
  50
  51 FilterTurbulence::~FilterTurbulence()
  52 {
  53     if (pix) {
  54         nr_pixblock_release(pix);
  55         delete pix;
  56     }
  57 }
  58
  59 void FilterTurbulence::set_baseFrequency(int axis, double freq){
  60     if (axis==0) XbaseFrequency=freq;
  61     if (axis==1) YbaseFrequency=freq;
  62 }
  63
  64 void FilterTurbulence::set_numOctaves(int num){
  65     numOctaves=num;
  66 }
  67
  68 void FilterTurbulence::set_seed(double s){
  69     seed=s;
  70 }
  71
  72 void FilterTurbulence::set_stitchTiles(bool st){
  73     stitchTiles=st;
  74 }
  75
  76 void FilterTurbulence::set_type(FilterTurbulenceType t){
  77     type=t;
  78 }
  79
  80 void FilterTurbulence::set_updated(bool u){
  81     updated=u;
  82 }
  83
  84 void FilterTurbulence::render_area(NRPixBlock *pix, NR::IRect &full_area, FilterUnits const &units) {
  85     const int bbox_x0 = full_area.min()[NR::X];
  86     const int bbox_y0 = full_area.min()[NR::Y];
  87     const int bbox_x1 = full_area.max()[NR::X];
  88     const int bbox_y1 = full_area.max()[NR::Y];
  89
  90     Geom::Matrix unit_trans = units.get_matrix_primitiveunits2pb().inverse();
  91
  92     double point[2];
  93
  94     unsigned char *pb = NR_PIXBLOCK_PX(pix);
  95
  96     if (type==TURBULENCE_TURBULENCE){
  97         for (int y = std::max(bbox_y0, pix->area.y0); y < std::min(bbox_y1, pix->area.y1); y++){
  98             int out_line = (y - pix->area.y0) * pix->rs;
  99             point[1] = y * unit_trans[3] + unit_trans[5];
 100             for (int x = std::max(bbox_x0, pix->area.x0); x < std::min(bbox_x1, pix->area.x1); x++){
 101                 int out_pos = out_line + 4 * (x - pix->area.x0);
 102                 point[0] = x * unit_trans[0] + unit_trans[4];
 103                 pb[out_pos] = CLAMP_D_TO_U8( turbulence(0,point)*255 ); // CLAMP includes rounding!
 104                 pb[out_pos + 1] = CLAMP_D_TO_U8( turbulence(1,point)*255 );
 105                 pb[out_pos + 2] = CLAMP_D_TO_U8( turbulence(2,point)*255 );
 106                 pb[out_pos + 3] = CLAMP_D_TO_U8( turbulence(3,point)*255 );
 107             }
 108         }
 109     } else {
 110         for (int y = std::max(bbox_y0, pix->area.y0); y < std::min(bbox_y1, pix->area.y1); y++){
 111             int out_line = (y - pix->area.y0) * pix->rs;
 112             point[1] = y * unit_trans[3] + unit_trans[5];
 113             for (int x = std::max(bbox_x0, pix->area.x0); x < std::min(bbox_x1, pix->area.x1); x++){
 114                 int out_pos = out_line + 4 * (x - pix->area.x0);
 115                 point[0] = x * unit_trans[0] + unit_trans[4];
 116                 pb[out_pos] = CLAMP_D_TO_U8( ((turbulence(0,point)*255) +255)/2 );
 117                 pb[out_pos + 1] = CLAMP_D_TO_U8( ((turbulence(1,point)*255)+255)/2 );
 118                 pb[out_pos + 2] = CLAMP_D_TO_U8( ((turbulence(2,point)*255) +255)/2 );
 119                 pb[out_pos + 3] = CLAMP_D_TO_U8( ((turbulence(3,point)*255) +255)/2 );
 120             }
 121         }
 122     }
 123
 124     pix->empty = FALSE;
 125 }
 126
 127 void FilterTurbulence::update_pixbuffer(NR::IRect &area, FilterUnits const &units) {
 128     int bbox_x0 = area.min()[NR::X];
 129     int bbox_y0 = area.min()[NR::Y];
 130     int bbox_x1 = area.max()[NR::X];
 131     int bbox_y1 = area.max()[NR::Y];
 132
 133     TurbulenceInit((long)seed);
 134
 135     if (!pix){
 136         pix = new NRPixBlock;
 137         nr_pixblock_setup_fast(pix, NR_PIXBLOCK_MODE_R8G8B8A8N, bbox_x0, bbox_y0, bbox_x1, bbox_y1, true);
 138     }
 139     else if (bbox_x0 != pix->area.x0 || bbox_y0 != pix->area.y0 ||
 140         bbox_x1 != pix->area.x1 || bbox_y1 != pix->area.y1)
 141     {
 142         /* TODO: release-setup cycle not actually needed, if pixblock
 143          * width and height don't change */
 144         nr_pixblock_release(pix);
 145         nr_pixblock_setup_fast(pix, NR_PIXBLOCK_MODE_R8G8B8A8N, bbox_x0, bbox_y0, bbox_x1, bbox_y1, true);
 146     }
 147
 148     /* This limits pre-rendered turbulence to two megapixels. This is
 149      * arbitary limit and could be something other, too.
 150      * If bigger area is needed, visible area is rendered on demand. */
 151     if (!pix || (pix->size != NR_PIXBLOCK_SIZE_TINY && pix->data.px == NULL) ||
 152         ((bbox_x1 - bbox_x0) * (bbox_y1 - bbox_y0) > 2*1024*1024)) {
 153         pix_data = NULL;
 154         return;
 155     }
 156
 157     render_area(pix, area, units);
 158
 159     pix_data = NR_PIXBLOCK_PX(pix);
 160
 161     updated=true;
 162     updated_area = area;
 163 }
 164
 165 int FilterTurbulence::render(FilterSlot &slot, FilterUnits const &units) {
 166     NR::IRect area = units.get_pixblock_filterarea_paraller();
 167     // TODO: could be faster - updated_area only has to be same size as area
 168     if (!updated || updated_area != area) update_pixbuffer(area, units);
 169
 170     NRPixBlock *in = slot.get(_input);
 171     if (!in) {
 172         g_warning("Missing source image for feTurbulence (in=%d)", _input);
 173         return 1;
 174     }
 175
 176     NRPixBlock *out = new NRPixBlock;
 177     int x0 = in->area.x0, y0 = in->area.y0;
 178     int x1 = in->area.x1, y1 = in->area.y1;
 179     nr_pixblock_setup_fast(out, NR_PIXBLOCK_MODE_R8G8B8A8N, x0, y0, x1, y1, true);
 180
 181     if (pix_data) {
 182         /* If pre-rendered output of whole filter area exists, just copy it. */
 183         nr_blit_pixblock_pixblock(out, pix);
 184     } else {
 185         /* No pre-rendered output, render the required area here. */
 186         render_area(out, area, units);
 187     }
 188
 189     out->empty = FALSE;
 190     slot.set(_output, out);
 191     return 0;
 192 }
 193
 194 long FilterTurbulence::Turbulence_setup_seed(long lSeed)
 195 {
 196   if (lSeed <= 0) lSeed = -(lSeed % (RAND_m - 1)) + 1;
 197   if (lSeed > RAND_m - 1) lSeed = RAND_m - 1;
 198   return lSeed;
 199 }
 200
 201 long FilterTurbulence::TurbulenceRandom(long lSeed)
 202 {
 203   long result;
 204   result = RAND_a * (lSeed % RAND_q) - RAND_r * (lSeed / RAND_q);
 205   if (result <= 0) result += RAND_m;
 206   return result;
 207 }
 208
 209 void FilterTurbulence::TurbulenceInit(long lSeed)
 210 {
 211   double s;
 212   int i, j, k;
 213   lSeed = Turbulence_setup_seed(lSeed);
 214   for(k = 0; k < 4; k++)
 215   {
 216     for(i = 0; i < BSize; i++)
 217     {
 218       uLatticeSelector[i] = i;
 219       for (j = 0; j < 2; j++)
 220         fGradient[k][i][j] = (double)(((lSeed = TurbulenceRandom(lSeed)) % (BSize + BSize)) - BSize) / BSize;
 221       s = double(sqrt(fGradient[k][i][0] * fGradient[k][i][0] + fGradient[k][i][1] * fGradient[k][i][1]));
 222       fGradient[k][i][0] /= s;
 223       fGradient[k][i][1] /= s;
 224     }
 225   }
 226   while(--i)
 227   {
 228     k = uLatticeSelector[i];
 229     uLatticeSelector[i] = uLatticeSelector[j = (lSeed = TurbulenceRandom(lSeed)) % BSize];
 230     uLatticeSelector[j] = k;
 231   }
 232   for(i = 0; i < BSize + 2; i++)
 233   {
 234     uLatticeSelector[BSize + i] = uLatticeSelector[i];
 235     for(k = 0; k < 4; k++)
 236       for(j = 0; j < 2; j++)
 237         fGradient[k][BSize + i][j] = fGradient[k][i][j];
 238   }
 239 }
 240
 241 double FilterTurbulence::TurbulenceNoise2(int nColorChannel, double vec[2], StitchInfo *pStitchInfo)
 242 {
 243   int bx0, bx1, by0, by1, b00, b10, b01, b11;
 244   double rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
 245   int i, j;
 246   t = vec[0] + PerlinN;
 247   bx0 = (int)t;
 248   bx1 = bx0+1;
 249   rx0 = t - (int)t;
 250   rx1 = rx0 - 1.0f;
 251   t = vec[1] + PerlinN;
 252   by0 = (int)t;
 253   by1 = by0+1;
 254   ry0 = t - (int)t;
 255   ry1 = ry0 - 1.0f;
 256   // If stitching, adjust lattice points accordingly.
 257   if(pStitchInfo != NULL)
 258   {
 259     if(bx0 >= pStitchInfo->nWrapX)
 260       bx0 -= pStitchInfo->nWidth;
 261     if(bx1 >= pStitchInfo->nWrapX)
 262       bx1 -= pStitchInfo->nWidth;
 263     if(by0 >= pStitchInfo->nWrapY)
 264       by0 -= pStitchInfo->nHeight;
 265     if(by1 >= pStitchInfo->nWrapY)
 266       by1 -= pStitchInfo->nHeight;
 267   }
 268   bx0 &= BM;
 269   bx1 &= BM;
 270   by0 &= BM;
 271   by1 &= BM;
 272   i = uLatticeSelector[bx0];
 273   j = uLatticeSelector[bx1];
 274   b00 = uLatticeSelector[i + by0];
 275   b10 = uLatticeSelector[j + by0];
 276   b01 = uLatticeSelector[i + by1];
 277   b11 = uLatticeSelector[j + by1];
 278   sx = double(s_curve(rx0));
 279   sy = double(s_curve(ry0));
 280   q = fGradient[nColorChannel][b00]; u = rx0 * q[0] + ry0 * q[1];
 281   q = fGradient[nColorChannel][b10]; v = rx1 * q[0] + ry0 * q[1];
 282   a = turb_lerp(sx, u, v);
 283   q = fGradient[nColorChannel][b01]; u = rx0 * q[0] + ry1 * q[1];
 284   q = fGradient[nColorChannel][b11]; v = rx1 * q[0] + ry1 * q[1];
 285   b = turb_lerp(sx, u, v);
 286   return turb_lerp(sy, a, b);
 287 }
 288
 289 double FilterTurbulence::turbulence(int nColorChannel, double *point)
 290 {
 291   StitchInfo stitch;
 292   StitchInfo *pStitchInfo = NULL; // Not stitching when NULL.
 293   // Adjust the base frequencies if necessary for stitching.
 294   if(stitchTiles)
 295   {
 296     // When stitching tiled turbulence, the frequencies must be adjusted
 297     // so that the tile borders will be continuous.
 298     if(XbaseFrequency != 0.0)
 299     {
 300       double fLoFreq = double(floor(fTileWidth * XbaseFrequency)) / fTileWidth;
 301       double fHiFreq = double(ceil(fTileWidth * XbaseFrequency)) / fTileWidth;
 302       if(XbaseFrequency / fLoFreq < fHiFreq / XbaseFrequency)
 303         XbaseFrequency = fLoFreq;
 304       else
 305         XbaseFrequency = fHiFreq;
 306     }
 307     if(YbaseFrequency != 0.0)
 308     {
 309       double fLoFreq = double(floor(fTileHeight * YbaseFrequency)) / fTileHeight;
 310       double fHiFreq = double(ceil(fTileHeight * YbaseFrequency)) / fTileHeight;
 311       if(YbaseFrequency / fLoFreq < fHiFreq / YbaseFrequency)
 312         YbaseFrequency = fLoFreq;
 313       else
 314         YbaseFrequency = fHiFreq;
 315     }
 316     // Set up TurbulenceInitial stitch values.
 317     pStitchInfo = &stitch;
 318     stitch.nWidth = int(fTileWidth * XbaseFrequency + 0.5f);
 319     stitch.nWrapX = int(fTileX * XbaseFrequency + PerlinN + stitch.nWidth);
 320     stitch.nHeight = int(fTileHeight * YbaseFrequency + 0.5f);
 321     stitch.nWrapY = int(fTileY * YbaseFrequency + PerlinN + stitch.nHeight);
 322   }
 323   double fSum = 0.0f;
 324   double vec[2];
 325   vec[0] = point[0] * XbaseFrequency;
 326   vec[1] = point[1] * YbaseFrequency;
 327   double ratio = 1;
 328   for(int nOctave = 0; nOctave < numOctaves; nOctave++)
 329   {
 330     if(type==TURBULENCE_FRACTALNOISE)
 331       fSum += double(TurbulenceNoise2(nColorChannel, vec, pStitchInfo) / ratio);
 332     else
 333       fSum += double(fabs(TurbulenceNoise2(nColorChannel, vec, pStitchInfo)) / ratio);
 334     vec[0] *= 2;
 335     vec[1] *= 2;
 336     ratio *= 2;
 337     if(pStitchInfo != NULL)
 338     {
 339       // Update stitch values. Subtracting PerlinN before the multiplication and
 340       // adding it afterward simplifies to subtracting it once.
 341       stitch.nWidth *= 2;
 342       stitch.nWrapX = 2 * stitch.nWrapX - PerlinN;
 343       stitch.nHeight *= 2;
 344       stitch.nWrapY = 2 * stitch.nWrapY - PerlinN;
 345     }
 346   }
 347   return fSum;
 348 }
 349
 350 FilterTraits FilterTurbulence::get_input_traits() {
 351     return TRAIT_PARALLER;
 352 }
 353
 354 } /* namespace Filters */
 355 } /* namespace Inkscape */
 356
 357 /*
 358   Local Variables:
 359   mode:c++
 360   c-file-style:"stroustrup"
 361   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 362   indent-tabs-mode:nil
 363   fill-column:99
 364   End:
 365 */
 366 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :