source/blender/compositor/intern/COM_MemoryBuffer.cpp

   1 /*
   2  * Copyright 2011, Blender Foundation.
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License
   6  * as published by the Free Software Foundation; either version 2
   7  * of the License, or (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program; if not, write to the Free Software Foundation,
  16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  17  *
  18  * Contributor:
  19  *              Jeroen Bakker
  20  *              Monique Dewanchand
  21  */
  22
  23 #include "COM_MemoryBuffer.h"
  24 #include "MEM_guardedalloc.h"
  25 //#include "BKE_global.h"
  26
  27 unsigned int MemoryBuffer::determineBufferSize()
  28 {
  29         return getWidth() * getHeight();
  30 }
  31
  32 int MemoryBuffer::getWidth() const
  33 {
  34         return this->m_rect.xmax - this->m_rect.xmin;
  35 }
  36 int MemoryBuffer::getHeight() const
  37 {
  38         return this->m_rect.ymax - this->m_rect.ymin;
  39 }
  40
  41 MemoryBuffer::MemoryBuffer(MemoryProxy *memoryProxy, unsigned int chunkNumber, rcti *rect)
  42 {
  43         BLI_rcti_init(&this->m_rect, rect->xmin, rect->xmax, rect->ymin, rect->ymax);
  44         this->m_memoryProxy = memoryProxy;
  45         this->m_chunkNumber = chunkNumber;
  46         this->m_buffer = (float *)MEM_mallocN(sizeof(float) * determineBufferSize() * COM_NUMBER_OF_CHANNELS, "COM_MemoryBuffer");
  47         this->m_state = COM_MB_ALLOCATED;
  48         this->m_datatype = COM_DT_COLOR;
  49         this->m_chunkWidth = this->m_rect.xmax - this->m_rect.xmin;
  50 }
  51
  52 MemoryBuffer::MemoryBuffer(MemoryProxy *memoryProxy, rcti *rect)
  53 {
  54         BLI_rcti_init(&this->m_rect, rect->xmin, rect->xmax, rect->ymin, rect->ymax);
  55         this->m_memoryProxy = memoryProxy;
  56         this->m_chunkNumber = -1;
  57         this->m_buffer = (float *)MEM_mallocN(sizeof(float) * determineBufferSize() * COM_NUMBER_OF_CHANNELS, "COM_MemoryBuffer");
  58         this->m_state = COM_MB_TEMPORARILY;
  59         this->m_datatype = COM_DT_COLOR;
  60         this->m_chunkWidth = this->m_rect.xmax - this->m_rect.xmin;
  61 }
  62 MemoryBuffer *MemoryBuffer::duplicate()
  63 {
  64         MemoryBuffer *result = new MemoryBuffer(this->m_memoryProxy, &this->m_rect);
  65         memcpy(result->m_buffer, this->m_buffer, this->determineBufferSize() * COM_NUMBER_OF_CHANNELS * sizeof(float));
  66         return result;
  67 }
  68 void MemoryBuffer::clear()
  69 {
  70         memset(this->m_buffer, 0, this->determineBufferSize() * COM_NUMBER_OF_CHANNELS * sizeof(float));
  71 }
  72
  73 float *MemoryBuffer::convertToValueBuffer()
  74 {
  75         const unsigned int size = this->determineBufferSize();
  76         unsigned int i;
  77
  78         float *result = (float *)MEM_mallocN(sizeof(float) * size, __func__);
  79
  80         const float *fp_src = this->m_buffer;
  81         float       *fp_dst = result;
  82
  83         for (i = 0; i < size; i++, fp_dst++, fp_src += COM_NUMBER_OF_CHANNELS) {
  84                 *fp_dst = *fp_src;
  85         }
  86
  87         return result;
  88 }
  89
  90 float MemoryBuffer::getMaximumValue()
  91 {
  92         float result = this->m_buffer[0];
  93         const unsigned int size = this->determineBufferSize();
  94         unsigned int i;
  95
  96         const float *fp_src = this->m_buffer;
  97
  98         for (i = 0; i < size; i++, fp_src += COM_NUMBER_OF_CHANNELS) {
  99                 float value = *fp_src;
 100                 if (value > result) {
 101                         result = value;
 102                 }
 103         }
 104
 105         return result;
 106 }
 107
 108 float MemoryBuffer::getMaximumValue(rcti *rect)
 109 {
 110         rcti rect_clamp;
 111
 112         /* first clamp the rect by the bounds or we get un-initialized values */
 113         BLI_rcti_isect(rect, &this->m_rect, &rect_clamp);
 114
 115         if (!BLI_rcti_is_empty(&rect_clamp)) {
 116                 MemoryBuffer *temp = new MemoryBuffer(NULL, &rect_clamp);
 117                 temp->copyContentFrom(this);
 118                 float result = temp->getMaximumValue();
 119                 delete temp;
 120                 return result;
 121         }
 122         else {
 123                 BLI_assert(0);
 124                 return 0.0f;
 125         }
 126 }
 127
 128 MemoryBuffer::~MemoryBuffer()
 129 {
 130         if (this->m_buffer) {
 131                 MEM_freeN(this->m_buffer);
 132                 this->m_buffer = NULL;
 133         }
 134 }
 135
 136 void MemoryBuffer::copyContentFrom(MemoryBuffer *otherBuffer)
 137 {
 138         if (!otherBuffer) {
 139                 BLI_assert(0);
 140                 return;
 141         }
 142         unsigned int otherY;
 143         unsigned int minX = max(this->m_rect.xmin, otherBuffer->m_rect.xmin);
 144         unsigned int maxX = min(this->m_rect.xmax, otherBuffer->m_rect.xmax);
 145         unsigned int minY = max(this->m_rect.ymin, otherBuffer->m_rect.ymin);
 146         unsigned int maxY = min(this->m_rect.ymax, otherBuffer->m_rect.ymax);
 147         int offset;
 148         int otherOffset;
 149
 150
 151         for (otherY = minY; otherY < maxY; otherY++) {
 152                 otherOffset = ((otherY - otherBuffer->m_rect.ymin) * otherBuffer->m_chunkWidth + minX - otherBuffer->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
 153                 offset = ((otherY - this->m_rect.ymin) * this->m_chunkWidth + minX - this->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
 154                 memcpy(&this->m_buffer[offset], &otherBuffer->m_buffer[otherOffset], (maxX - minX) * COM_NUMBER_OF_CHANNELS * sizeof(float));
 155         }
 156 }
 157
 158 void MemoryBuffer::writePixel(int x, int y, const float color[4])
 159 {
 160         if (x >= this->m_rect.xmin && x < this->m_rect.xmax &&
 161             y >= this->m_rect.ymin && y < this->m_rect.ymax)
 162         {
 163                 const int offset = (this->m_chunkWidth * (y-this->m_rect.ymin) + x-this->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
 164                 copy_v4_v4(&this->m_buffer[offset], color);
 165         }
 166 }
 167
 168 void MemoryBuffer::addPixel(int x, int y, const float color[4])
 169 {
 170         if (x >= this->m_rect.xmin && x < this->m_rect.xmax &&
 171             y >= this->m_rect.ymin && y < this->m_rect.ymax)
 172         {
 173                 const int offset = (this->m_chunkWidth * (y-this->m_rect.ymin) + x-this->m_rect.xmin) * COM_NUMBER_OF_CHANNELS;
 174                 add_v4_v4(&this->m_buffer[offset], color);
 175         }
 176 }
 177
 178
 179 // table of (exp(ar) - exp(a)) / (1 - exp(a)) for r in range [0, 1] and a = -2
 180 // used instead of actual gaussian, otherwise at high texture magnifications circular artifacts are visible
 181 #define EWA_MAXIDX 255
 182 static float EWA_WTS[EWA_MAXIDX + 1] = {
 183         1.f, 0.990965f, 0.982f, 0.973105f, 0.96428f, 0.955524f, 0.946836f, 0.938216f, 0.929664f,
 184         0.921178f, 0.912759f, 0.904405f, 0.896117f, 0.887893f, 0.879734f, 0.871638f, 0.863605f,
 185         0.855636f, 0.847728f, 0.839883f, 0.832098f, 0.824375f, 0.816712f, 0.809108f, 0.801564f,
 186         0.794079f, 0.786653f, 0.779284f, 0.771974f, 0.76472f, 0.757523f, 0.750382f, 0.743297f,
 187         0.736267f, 0.729292f, 0.722372f, 0.715505f, 0.708693f, 0.701933f, 0.695227f, 0.688572f,
 188         0.68197f, 0.67542f, 0.66892f, 0.662471f, 0.656073f, 0.649725f, 0.643426f, 0.637176f,
 189         0.630976f, 0.624824f, 0.618719f, 0.612663f, 0.606654f, 0.600691f, 0.594776f, 0.588906f,
 190         0.583083f, 0.577305f, 0.571572f, 0.565883f, 0.56024f, 0.55464f, 0.549084f, 0.543572f,
 191         0.538102f, 0.532676f, 0.527291f, 0.521949f, 0.516649f, 0.511389f, 0.506171f, 0.500994f,
 192         0.495857f, 0.490761f, 0.485704f, 0.480687f, 0.475709f, 0.470769f, 0.465869f, 0.461006f,
 193         0.456182f, 0.451395f, 0.446646f, 0.441934f, 0.437258f, 0.432619f, 0.428017f, 0.42345f,
 194         0.418919f, 0.414424f, 0.409963f, 0.405538f, 0.401147f, 0.39679f, 0.392467f, 0.388178f,
 195         0.383923f, 0.379701f, 0.375511f, 0.371355f, 0.367231f, 0.363139f, 0.359079f, 0.355051f,
 196         0.351055f, 0.347089f, 0.343155f, 0.339251f, 0.335378f, 0.331535f, 0.327722f, 0.323939f,
 197         0.320186f, 0.316461f, 0.312766f, 0.3091f, 0.305462f, 0.301853f, 0.298272f, 0.294719f,
 198         0.291194f, 0.287696f, 0.284226f, 0.280782f, 0.277366f, 0.273976f, 0.270613f, 0.267276f,
 199         0.263965f, 0.26068f, 0.257421f, 0.254187f, 0.250979f, 0.247795f, 0.244636f, 0.241502f,
 200         0.238393f, 0.235308f, 0.232246f, 0.229209f, 0.226196f, 0.223206f, 0.220239f, 0.217296f,
 201         0.214375f, 0.211478f, 0.208603f, 0.20575f, 0.20292f, 0.200112f, 0.197326f, 0.194562f,
 202         0.191819f, 0.189097f, 0.186397f, 0.183718f, 0.18106f, 0.178423f, 0.175806f, 0.17321f,
 203         0.170634f, 0.168078f, 0.165542f, 0.163026f, 0.16053f, 0.158053f, 0.155595f, 0.153157f,
 204         0.150738f, 0.148337f, 0.145955f, 0.143592f, 0.141248f, 0.138921f, 0.136613f, 0.134323f,
 205         0.132051f, 0.129797f, 0.12756f, 0.125341f, 0.123139f, 0.120954f, 0.118786f, 0.116635f,
 206         0.114501f, 0.112384f, 0.110283f, 0.108199f, 0.106131f, 0.104079f, 0.102043f, 0.100023f,
 207         0.0980186f, 0.09603f, 0.094057f, 0.0920994f, 0.0901571f, 0.08823f, 0.0863179f, 0.0844208f,
 208         0.0825384f, 0.0806708f, 0.0788178f, 0.0769792f, 0.0751551f, 0.0733451f, 0.0715493f, 0.0697676f,
 209         0.0679997f, 0.0662457f, 0.0645054f, 0.0627786f, 0.0610654f, 0.0593655f, 0.0576789f, 0.0560055f,
 210         0.0543452f, 0.0526979f, 0.0510634f, 0.0494416f, 0.0478326f, 0.0462361f, 0.0446521f, 0.0430805f,
 211         0.0415211f, 0.039974f, 0.0384389f, 0.0369158f, 0.0354046f, 0.0339052f, 0.0324175f, 0.0309415f,
 212         0.029477f, 0.0280239f, 0.0265822f, 0.0251517f, 0.0237324f, 0.0223242f, 0.020927f, 0.0195408f,
 213         0.0181653f, 0.0168006f, 0.0154466f, 0.0141031f, 0.0127701f, 0.0114476f, 0.0101354f, 0.00883339f,
 214         0.00754159f, 0.00625989f, 0.00498819f, 0.00372644f, 0.00247454f, 0.00123242f, 0.f
 215 };
 216
 217 static void radangle2imp(float a2, float b2, float th, float *A, float *B, float *C, float *F)
 218 {
 219         float ct2 = cosf(th);
 220         const float st2 = 1.f - ct2 * ct2;    // <- sin(th)^2
 221         ct2 *= ct2;
 222         *A = a2 * st2 + b2 * ct2;
 223         *B = (b2 - a2) * sinf(2.f * th);
 224         *C = a2 * ct2 + b2 * st2;
 225         *F = a2 * b2;
 226 }
 227
 228 // all tests here are done to make sure possible overflows are hopefully minimized
 229 static void imp2radangle(float A, float B, float C, float F, float *a, float *b, float *th, float *ecc)
 230 {
 231         if (F <= 1e-5f) {   // use arbitrary major radius, zero minor, infinite eccentricity
 232                 *a = sqrtf(A > C ? A : C);
 233                 *b = 0.f;
 234                 *ecc = 1e10f;
 235                 *th = 0.5f * (atan2f(B, A - C) + (float)M_PI);
 236         }
 237         else {
 238                 const float AmC = A - C, ApC = A + C, F2 = F * 2.f;
 239                 const float r = sqrtf(AmC * AmC + B * B);
 240                 float d = ApC - r;
 241                 *a = (d <= 0.f) ? sqrtf(A > C ? A : C) : sqrtf(F2 / d);
 242                 d = ApC + r;
 243                 if (d <= 0.f) {
 244                         *b = 0.f;
 245                         *ecc = 1e10f;
 246                 }
 247                 else {
 248                         *b = sqrtf(F2 / d);
 249                         *ecc = *a / *b;
 250                 }
 251                 /* incr theta by 0.5 * pi (angle of major axis) */
 252                 *th = 0.5f * (atan2f(B, AmC) + (float)M_PI);
 253         }
 254 }
 255
 256 float clipuv(float x, float limit)
 257 {
 258         x = (x < 0) ? 0 : ((x >= limit) ? (limit - 1) : x);
 259         return x;
 260 }
 261
 262 /**
 263  * \note \a sampler at the moment is either 'COM_PS_NEAREST' or not, other values won't matter.
 264  */
 265 void MemoryBuffer::readEWA(float result[4], float fx, float fy, float dx, float dy, PixelSampler sampler)
 266 {
 267         const int width = this->getWidth(), height = this->getHeight();
 268
 269         // scaling dxt/dyt by full resolution can cause overflow because of huge A/B/C and esp. F values,
 270         // scaling by aspect ratio alone does the opposite, so try something in between instead...
 271         const float ff2 = width, ff = sqrtf(ff2), q = height / ff;
 272         const float Ux = dx * ff, Vx = dx * q, Uy = dy * ff, Vy = dy * q;
 273         float A = Vx * Vx + Vy * Vy;
 274         float B = -2.f * (Ux * Vx + Uy * Vy);
 275         float C = Ux * Ux + Uy * Uy;
 276         float F = A * C - B * B * 0.25f;
 277         float a, b, th, ecc, a2, b2, ue, ve, U0, V0, DDQ, U, ac1, ac2, BU, d;
 278         int u, v, u1, u2, v1, v2;
 279         // The so-called 'high' quality ewa method simply adds a constant of 1 to both A & C,
 280         // so the ellipse always covers at least some texels. But since the filter is now always larger,
 281         // it also means that everywhere else it's also more blurry then ideally should be the case.
 282         // So instead here the ellipse radii are modified instead whenever either is too low.
 283         // Use a different radius based on interpolation switch, just enough to anti-alias when interpolation is off,
 284         // and slightly larger to make result a bit smoother than bilinear interpolation when interpolation is on
 285         // (minimum values: const float rmin = intpol ? 1.f : 0.5f;)
 286
 287         /* note: 0.765625f is too sharp, 1.0 will not blur with an exact pixel sample
 288          * useful to avoid blurring when there is no distortion */
 289 #if 0
 290         const float rmin = ((sampler != COM_PS_NEAREST) ? 1.5625f : 0.765625f) / ff2;
 291 #else
 292         const float rmin = ((sampler != COM_PS_NEAREST) ? 1.5625f : 1.0f     ) / ff2;
 293 #endif
 294         imp2radangle(A, B, C, F, &a, &b, &th, &ecc);
 295         if ((b2 = b * b) < rmin) {
 296                 if ((a2 = a * a) < rmin) {
 297                         B = 0.f;
 298                         A = C = rmin;
 299                         F = A * C;
 300                 }
 301                 else {
 302                         b2 = rmin;
 303                         radangle2imp(a2, b2, th, &A, &B, &C, &F);
 304                 }
 305         }
 306
 307         ue = ff * sqrtf(C);
 308         ve = ff * sqrtf(A);
 309         d = (float)(EWA_MAXIDX + 1) / (F * ff2);
 310         A *= d;
 311         B *= d;
 312         C *= d;
 313
 314         U0 = fx;
 315         V0 = fy;
 316         u1 = (int)(floorf(U0 - ue));
 317         u2 = (int)(ceilf(U0 + ue));
 318         v1 = (int)(floorf(V0 - ve));
 319         v2 = (int)(ceilf(V0 + ve));
 320         U0 -= 0.5f;
 321         V0 -= 0.5f;
 322         DDQ = 2.f * A;
 323         U = u1 - U0;
 324         ac1 = A * (2.f * U + 1.f);
 325         ac2 = A * U * U;
 326         BU = B * U;
 327
 328         d = result[0] = result[1] = result[2] = result[3] = 0.f;
 329         for (v = v1; v <= v2; ++v) {
 330                 const float V = v - V0;
 331                 float DQ = ac1 + B * V;
 332                 float Q = (C * V + BU) * V + ac2;
 333                 for (u = u1; u <= u2; ++u) {
 334                         if (Q < (float)(EWA_MAXIDX + 1)) {
 335                                 float tc[4];
 336                                 const float wt = EWA_WTS[(Q < 0.f) ? 0 : (unsigned int)Q];
 337                                 read(tc, clipuv(u, width), clipuv(v, height));
 338                                 madd_v3_v3fl(result, tc, wt);
 339                                 result[3] += result[3] ? tc[3] * wt : 0.f;
 340                                 d += wt;
 341                         }
 342                         Q += DQ;
 343                         DQ += DDQ;
 344                 }
 345         }
 346
 347         // d should hopefully never be zero anymore
 348         d = 1.f / d;
 349         result[0] *= d;
 350         result[1] *= d;
 351         result[2] *= d;
 352         // clipping can be ignored if alpha used, texr->ta already includes filtered edge
 353         result[3] = result[3] ? result[3] * d : 1.f;
 354 }