1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
2
3 /*
4 Sonic Visualiser
5 An audio file viewer and annotation editor.
6 Centre for Digital Music, Queen Mary, University of London.
7 This file copyright 2006 Chris Cannam and QMUL.
8
9 This program is free software; you can redistribute it and/or
10 modify it under the terms of the GNU General Public License as
11 published by the Free Software Foundation; either version 2 of the
12 License, or (at your option) any later version. See the file
13 COPYING included with this distribution for more information.
14 */
15
16 #include "AudioTimeStretcher.h"
17
18 #include <QtGlobal>
19
20 #include <fstream>
21 #include <cstring>
22
23 namespace Rosegarden
24 {
25
26 // (Unused)
27 // static double mod(double x, double y) { return x - (y * floor(x / y)); }
modf(float x,float y)28 static float modf(float x, float y) { return x - (y * floorf(x / y)); }
29
30 // (Unused)
31 // static double princarg(double a) { return mod(a + M_PI, -2 * M_PI) + M_PI; }
princargf(float a)32 static float princargf(float a) { return modf(a + M_PI, -2 * M_PI) + M_PI; }
33
34
35 //#define DEBUG_AUDIO_TIME_STRETCHER 1
36
AudioTimeStretcher(size_t sampleRate,size_t channels,float ratio,bool sharpen,size_t maxOutputBlockSize)37 AudioTimeStretcher::AudioTimeStretcher(size_t sampleRate,
38 size_t channels,
39 float ratio,
40 bool sharpen,
41 size_t maxOutputBlockSize) :
42 m_sampleRate(sampleRate),
43 m_channels(channels),
44 m_maxOutputBlockSize(maxOutputBlockSize),
45 m_ratio(ratio),
46 m_sharpen(sharpen),
47 m_totalCount(0),
48 m_transientCount(0),
49 m_n2sum(0),
50 m_n2total(0),
51 m_adjustCount(50)
52 {
53 pthread_mutex_t initialisingMutex = PTHREAD_MUTEX_INITIALIZER;
54 memcpy(&m_mutex, &initialisingMutex, sizeof(pthread_mutex_t));
55
56 initialise();
57 }
58
~AudioTimeStretcher()59 AudioTimeStretcher::~AudioTimeStretcher()
60 {
61 std::cerr << "AudioTimeStretcher::~AudioTimeStretcher" << std::endl;
62
63 std::cerr << "AudioTimeStretcher::~AudioTimeStretcher: actual ratio = " << (m_totalCount > 0 ? (float (m_n2total) / float(m_totalCount * m_n1)) : 1.f) << ", ideal = " << m_ratio << ", nominal = " << getRatio() << ")" << std::endl;
64
65 cleanup();
66
67 pthread_mutex_destroy(&m_mutex);
68 }
69
70 void
initialise()71 AudioTimeStretcher::initialise()
72 {
73 std::cerr << "AudioTimeStretcher::initialise" << std::endl;
74
75 calculateParameters();
76
77 m_analysisWindow = new SampleWindow<float>(SampleWindow<float>::Hanning, m_wlen);
78 m_synthesisWindow = new SampleWindow<float>(SampleWindow<float>::Hanning, m_wlen);
79
80 m_prevPhase = new float *[m_channels];
81 m_prevAdjustedPhase = new float *[m_channels];
82
83 m_prevTransientMag = (float *)fftwf_malloc(sizeof(float) * (m_wlen / 2 + 1));
84 m_prevTransientScore = 0;
85 m_prevTransient = false;
86
87 m_tempbuf = (float *)fftwf_malloc(sizeof(float) * m_wlen);
88
89 m_time = new float *[m_channels];
90 m_freq = new fftwf_complex *[m_channels];
91 m_plan = new fftwf_plan[m_channels];
92 m_iplan = new fftwf_plan[m_channels];
93
94 m_inbuf = new RingBuffer<float> *[m_channels];
95 m_outbuf = new RingBuffer<float> *[m_channels];
96 m_mashbuf = new float *[m_channels];
97
98 m_modulationbuf = (float *)fftwf_malloc(sizeof(float) * m_wlen);
99
100 for (size_t c = 0; c < m_channels; ++c) {
101
102 m_prevPhase[c] = (float *)fftwf_malloc(sizeof(float) * (m_wlen / 2 + 1));
103 m_prevAdjustedPhase[c] = (float *)fftwf_malloc(sizeof(float) * (m_wlen / 2 + 1));
104
105 m_time[c] = (float *)fftwf_malloc(sizeof(float) * m_wlen);
106 m_freq[c] = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) *
107 (m_wlen / 2 + 1));
108
109 m_plan[c] = fftwf_plan_dft_r2c_1d(m_wlen, m_time[c], m_freq[c], FFTW_ESTIMATE);
110 m_iplan[c] = fftwf_plan_dft_c2r_1d(m_wlen, m_freq[c], m_time[c], FFTW_ESTIMATE);
111
112 m_outbuf[c] = new RingBuffer<float>
113 ((m_maxOutputBlockSize + m_wlen) * 2);
114 m_inbuf[c] = new RingBuffer<float>
115 (lrintf(m_outbuf[c]->getSize() / m_ratio) + m_wlen);
116
117 std::cerr << "making inbuf size " << m_inbuf[c]->getSize() << " (outbuf size is " << m_outbuf[c]->getSize() << ", ratio " << m_ratio << ")" << std::endl;
118
119
120 m_mashbuf[c] = (float *)fftwf_malloc(sizeof(float) * m_wlen);
121
122 for (size_t i = 0; i < m_wlen; ++i) {
123 m_mashbuf[c][i] = 0.0;
124 }
125
126 for (size_t i = 0; i <= m_wlen/2; ++i) {
127 m_prevPhase[c][i] = 0.0;
128 m_prevAdjustedPhase[c][i] = 0.0;
129 }
130 }
131
132 for (size_t i = 0; i < m_wlen; ++i) {
133 m_modulationbuf[i] = 0.0;
134 }
135
136 for (size_t i = 0; i <= m_wlen/2; ++i) {
137 m_prevTransientMag[i] = 0.0;
138 }
139 }
140
141 void
calculateParameters()142 AudioTimeStretcher::calculateParameters()
143 {
144 std::cerr << "AudioTimeStretcher::calculateParameters" << std::endl;
145
146 m_wlen = 1024;
147
148 //!!! In transient sharpening mode, we need to pick the window
149 //length so as to be more or less fixed in audio duration (i.e. we
150 //need to exploit the sample rate)
151
152 //!!! have to work out the relationship between wlen and transient
153 //threshold
154
155 if (m_ratio < 1) {
156 if (m_ratio < 0.4) {
157 m_n1 = 1024;
158 m_wlen = 2048;
159 } else if (m_ratio < 0.8) {
160 m_n1 = 512;
161 } else {
162 m_n1 = 256;
163 }
164 if (shouldSharpen()) {
165 m_wlen = 2048;
166 }
167 m_n2 = lrintf(m_n1 * m_ratio);
168 } else {
169 if (m_ratio > 2) {
170 m_n2 = 512;
171 m_wlen = 4096;
172 } else if (m_ratio > 1.6) {
173 m_n2 = 384;
174 m_wlen = 2048;
175 } else {
176 m_n2 = 256;
177 }
178 if (shouldSharpen()) {
179 if (m_wlen < 2048) m_wlen = 2048;
180 }
181 m_n1 = lrintf(m_n2 / m_ratio);
182 if (m_n1 == 0) {
183 m_n1 = 1;
184 m_n2 = m_ratio;
185 }
186 }
187
188 m_transientThreshold = lrintf(m_wlen / 4.5);
189
190 m_totalCount = 0;
191 m_transientCount = 0;
192 m_n2sum = 0;
193 m_n2total = 0;
194 m_n2list.clear();
195
196 std::cerr << "AudioTimeStretcher: channels = " << m_channels
197 << ", ratio = " << m_ratio
198 << ", n1 = " << m_n1 << ", n2 = " << m_n2 << ", wlen = "
199 << m_wlen << ", max = " << m_maxOutputBlockSize << std::endl;
200 // << ", outbuflen = " << m_outbuf[0]->getSize() << std::endl;
201 }
202
203 void
cleanup()204 AudioTimeStretcher::cleanup()
205 {
206 std::cerr << "AudioTimeStretcher::cleanup" << std::endl;
207
208 for (size_t c = 0; c < m_channels; ++c) {
209
210 fftwf_destroy_plan(m_plan[c]);
211 fftwf_destroy_plan(m_iplan[c]);
212
213 fftwf_free(m_time[c]);
214 fftwf_free(m_freq[c]);
215
216 fftwf_free(m_mashbuf[c]);
217 fftwf_free(m_prevPhase[c]);
218 fftwf_free(m_prevAdjustedPhase[c]);
219
220 delete m_inbuf[c];
221 delete m_outbuf[c];
222 }
223
224 fftwf_free(m_tempbuf);
225 fftwf_free(m_modulationbuf);
226 fftwf_free(m_prevTransientMag);
227
228 delete[] m_prevPhase;
229 delete[] m_prevAdjustedPhase;
230 delete[] m_inbuf;
231 delete[] m_outbuf;
232 delete[] m_mashbuf;
233 delete[] m_time;
234 delete[] m_freq;
235 delete[] m_plan;
236 delete[] m_iplan;
237
238 delete m_analysisWindow;
239 delete m_synthesisWindow;
240 }
241
242 void
setRatio(float ratio)243 AudioTimeStretcher::setRatio(float ratio)
244 {
245 pthread_mutex_lock(&m_mutex);
246
247 size_t formerWlen = m_wlen;
248 m_ratio = ratio;
249
250 std::cerr << "AudioTimeStretcher::setRatio: new ratio " << ratio
251 << std::endl;
252
253 calculateParameters();
254
255 if (m_wlen == formerWlen) {
256
257 // This is the only container whose size depends on m_ratio
258
259 RingBuffer<float> **newin = new RingBuffer<float> *[m_channels];
260
261 size_t formerSize = m_inbuf[0]->getSize();
262 size_t newSize = lrintf(m_outbuf[0]->getSize() / m_ratio) + m_wlen;
263
264 std::cerr << "resizing inbuf from " << formerSize << " to "
265 << newSize << " (outbuf size is " << m_outbuf[0]->getSize() << ", ratio " << m_ratio << ")" << std::endl;
266
267 if (formerSize != newSize) {
268
269 size_t ready = m_inbuf[0]->getReadSpace();
270
271 for (size_t c = 0; c < m_channels; ++c) {
272 newin[c] = new RingBuffer<float>(newSize);
273 }
274
275 if (ready > 0) {
276
277 size_t copy = std::min(ready, newSize);
278 float *tmp = new float[ready];
279
280 for (size_t c = 0; c < m_channels; ++c) {
281 m_inbuf[c]->read(tmp, ready);
282 newin[c]->write(tmp + ready - copy, copy);
283 }
284
285 delete[] tmp;
286 }
287
288 for (size_t c = 0; c < m_channels; ++c) {
289 delete m_inbuf[c];
290 }
291
292 delete[] m_inbuf;
293 m_inbuf = newin;
294 }
295
296 } else {
297
298 std::cerr << "wlen changed" << std::endl;
299 cleanup();
300 initialise();
301 }
302
303 pthread_mutex_unlock(&m_mutex);
304 }
305
306 size_t
getProcessingLatency() const307 AudioTimeStretcher::getProcessingLatency() const
308 {
309 return getWindowSize() - getInputIncrement();
310 }
311
312 size_t
getRequiredInputSamples() const313 AudioTimeStretcher::getRequiredInputSamples() const
314 {
315 size_t rv;
316 pthread_mutex_lock(&m_mutex);
317
318 if (m_inbuf[0]->getReadSpace() >= m_wlen) rv = 0;
319 else rv = m_wlen - m_inbuf[0]->getReadSpace();
320
321 pthread_mutex_unlock(&m_mutex);
322 return rv;
323 }
324
325 void
putInput(float ** input,size_t samples)326 AudioTimeStretcher::putInput(float **input, size_t samples)
327 {
328 pthread_mutex_lock(&m_mutex);
329
330 // We need to add samples from input to our internal buffer. When
331 // we have m_windowSize samples in the buffer, we can process it,
332 // move the samples back by m_n1 and write the output onto our
333 // internal output buffer. If we have (samples * ratio) samples
334 // in that, we can write m_n2 of them back to output and return
335 // (otherwise we have to write zeroes).
336
337 // When we process, we write m_wlen to our fixed output buffer
338 // (m_mashbuf). We then pull out the first m_n2 samples from that
339 // buffer, push them into the output ring buffer, and shift
340 // m_mashbuf left by that amount.
341
342 // The processing latency is then m_wlen - m_n2.
343
344 size_t consumed = 0;
345
346 while (consumed < samples) {
347
348 size_t writable = m_inbuf[0]->getWriteSpace();
349 writable = std::min(writable, samples - consumed);
350
351 if (writable == 0) {
352 #ifdef DEBUG_AUDIO_TIME_STRETCHER
353 std::cerr << "WARNING: AudioTimeStretcher::putInput: writable == 0 (inbuf has " << m_inbuf[0]->getReadSpace() << " samples available for reading, space for " << m_inbuf[0]->getWriteSpace() << " more)" << std::endl;
354 #endif
355 if (m_inbuf[0]->getReadSpace() < m_wlen ||
356 m_outbuf[0]->getWriteSpace() < m_n2) {
357 std::cerr << "WARNING: AudioTimeStretcher::putInput: Inbuf has " << m_inbuf[0]->getReadSpace() << ", outbuf has space for " << m_outbuf[0]->getWriteSpace() << " (n2 = " << m_n2 << ", wlen = " << m_wlen << "), won't be able to process" << std::endl;
358 break;
359 }
360 } else {
361
362 #ifdef DEBUG_AUDIO_TIME_STRETCHER
363 std::cerr << "writing " << writable << " from index " << consumed << " to inbuf, consumed will be " << consumed + writable << std::endl;
364 #endif
365
366 for (size_t c = 0; c < m_channels; ++c) {
367 m_inbuf[c]->write(input[c] + consumed, writable);
368 }
369 consumed += writable;
370 }
371
372 while (m_inbuf[0]->getReadSpace() >= m_wlen &&
373 m_outbuf[0]->getWriteSpace() >= m_n2) {
374
375 // We know we have at least m_wlen samples available
376 // in m_inbuf. We need to peek m_wlen of them for
377 // processing, and then read m_n1 to advance the read
378 // pointer.
379
380 for (size_t c = 0; c < m_channels; ++c) {
381
382 #ifndef NDEBUG
383 size_t got = m_inbuf[c]->peek(m_tempbuf, m_wlen);
384 Q_ASSERT(got == m_wlen);
385 #endif
386
387 analyseBlock(c, m_tempbuf);
388 }
389
390 bool transient = false;
391 if (shouldSharpen()) transient = isTransient();
392
393 size_t n2 = m_n2;
394
395 if (transient) {
396 n2 = m_n1;
397 }
398
399 ++m_totalCount;
400 if (transient) ++m_transientCount;
401
402 m_n2sum += n2;
403 m_n2total += n2;
404
405 if (m_totalCount > 50 && m_transientCount < m_totalCount) {
406
407 int fixed = m_transientCount * m_n1;
408
409 float idealTotal = m_totalCount * m_n1 * m_ratio;
410 float idealSquashy = idealTotal - fixed;
411
412 float squashyCount = m_totalCount - m_transientCount;
413
414 float fn2 = idealSquashy / squashyCount;
415
416 n2 = int(fn2);
417
418 float remainder = fn2 - n2;
419 if (drand48() < remainder) ++n2;
420
421 #ifdef DEBUG_AUDIO_TIME_STRETCHER
422 if (n2 != m_n2) {
423 std::cerr << m_n2 << " -> " << n2 << " (ideal = " << (idealSquashy / squashyCount) << ")" << std::endl;
424 }
425 #endif
426 }
427
428 for (size_t c = 0; c < m_channels; ++c) {
429
430 synthesiseBlock(c, m_mashbuf[c],
431 c == 0 ? m_modulationbuf : nullptr,
432 m_prevTransient ? m_n1 : m_n2);
433
434
435 #ifdef DEBUG_AUDIO_TIME_STRETCHER
436 std::cerr << "writing first " << m_n2 << " from mashbuf, skipping " << m_n1 << " on inbuf " << std::endl;
437 #endif
438 m_inbuf[c]->skip(m_n1);
439
440 for (size_t i = 0; i < n2; ++i) {
441 if (m_modulationbuf[i] > 0.f) {
442 m_mashbuf[c][i] /= m_modulationbuf[i];
443 }
444 }
445
446 m_outbuf[c]->write(m_mashbuf[c], n2);
447
448 for (size_t i = 0; i < m_wlen - n2; ++i) {
449 m_mashbuf[c][i] = m_mashbuf[c][i + n2];
450 }
451
452 for (size_t i = m_wlen - n2; i < m_wlen; ++i) {
453 m_mashbuf[c][i] = 0.0f;
454 }
455 }
456
457 m_prevTransient = transient;
458
459 for (size_t i = 0; i < m_wlen - n2; ++i) {
460 m_modulationbuf[i] = m_modulationbuf[i + n2];
461 }
462
463 for (size_t i = m_wlen - n2; i < m_wlen; ++i) {
464 m_modulationbuf[i] = 0.0f;
465 }
466
467 if (!transient) m_n2 = n2;
468 }
469
470
471 #ifdef DEBUG_AUDIO_TIME_STRETCHER
472 std::cerr << "loop ended: inbuf read space " << m_inbuf[0]->getReadSpace() << ", outbuf write space " << m_outbuf[0]->getWriteSpace() << std::endl;
473 #endif
474 }
475
476 #ifdef DEBUG_AUDIO_TIME_STRETCHER
477 std::cerr << "AudioTimeStretcher::putInput returning" << std::endl;
478 #endif
479
480 pthread_mutex_unlock(&m_mutex);
481
482 // std::cerr << "ratio: nominal: " << getRatio() << " actual: "
483 // << m_total2 << "/" << m_total1 << " = " << float(m_total2) / float(m_total1) << " ideal: " << m_ratio << std::endl;
484 }
485
486 size_t
getAvailableOutputSamples() const487 AudioTimeStretcher::getAvailableOutputSamples() const
488 {
489 pthread_mutex_lock(&m_mutex);
490
491 size_t rv = m_outbuf[0]->getReadSpace();
492
493 pthread_mutex_unlock(&m_mutex);
494 return rv;
495 }
496
497 void
getOutput(float ** output,size_t samples)498 AudioTimeStretcher::getOutput(float **output, size_t samples)
499 {
500 pthread_mutex_lock(&m_mutex);
501
502 if (m_outbuf[0]->getReadSpace() < samples) {
503 std::cerr << "WARNING: AudioTimeStretcher::getOutput: not enough data (yet?) (" << m_outbuf[0]->getReadSpace() << " < " << samples << ")" << std::endl;
504 size_t fill = samples - m_outbuf[0]->getReadSpace();
505 for (size_t c = 0; c < m_channels; ++c) {
506 for (size_t i = 0; i < fill; ++i) {
507 output[c][i] = 0.0;
508 }
509 m_outbuf[c]->read(output[c] + fill, m_outbuf[c]->getReadSpace());
510 }
511 } else {
512 #ifdef DEBUG_AUDIO_TIME_STRETCHER
513 std::cerr << "enough data - writing " << samples << " from outbuf" << std::endl;
514 #endif
515 for (size_t c = 0; c < m_channels; ++c) {
516 m_outbuf[c]->read(output[c], samples);
517 }
518 }
519
520 #ifdef DEBUG_AUDIO_TIME_STRETCHER
521 std::cerr << "AudioTimeStretcher::getOutput returning" << std::endl;
522 #endif
523
524 pthread_mutex_unlock(&m_mutex);
525 }
526
527 void
analyseBlock(size_t c,float * buf)528 AudioTimeStretcher::analyseBlock(size_t c, float *buf)
529 {
530 size_t i;
531
532 // buf contains m_wlen samples
533
534 #ifdef DEBUG_AUDIO_TIME_STRETCHER
535 std::cerr << "AudioTimeStretcher::analyseBlock (channel " << c << ")" << std::endl;
536 #endif
537
538 m_analysisWindow->cut(buf);
539
540 for (i = 0; i < m_wlen/2; ++i) {
541 float temp = buf[i];
542 buf[i] = buf[i + m_wlen/2];
543 buf[i + m_wlen/2] = temp;
544 }
545
546 for (i = 0; i < m_wlen; ++i) {
547 m_time[c][i] = buf[i];
548 }
549
550 fftwf_execute(m_plan[c]); // m_time -> m_freq
551 }
552
553 bool
isTransient()554 AudioTimeStretcher::isTransient()
555 {
556 int count = 0;
557
558 for (size_t i = 0; i <= m_wlen/2; ++i) {
559
560 float real = 0.f, imag = 0.f;
561
562 for (size_t c = 0; c < m_channels; ++c) {
563 real += m_freq[c][i][0];
564 imag += m_freq[c][i][1];
565 }
566
567 float sqrmag = (real * real + imag * imag);
568
569 if (m_prevTransientMag[i] > 0.f) {
570 float diff = 10.f * log10f(sqrmag / m_prevTransientMag[i]);
571 if (diff > 3.f) ++count;
572 }
573
574 m_prevTransientMag[i] = sqrmag;
575 }
576
577 bool isTransient = false;
578
579 // if (count > m_transientThreshold &&
580 // count > m_prevTransientScore * 1.2) {
581 if (count > m_prevTransientScore &&
582 count > m_transientThreshold &&
583 int(count) - m_prevTransientScore > int(m_wlen) / 20) {
584 isTransient = true;
585
586 #ifdef DEBUG_AUDIO_TIME_STRETCHER
587 std::cerr << "isTransient (count = " << count << ", prev = " << m_prevTransientScore << ", diff = " << count - m_prevTransientScore << ", ratio = " << (m_totalCount > 0 ? (float (m_n2sum) / float(m_totalCount * m_n1)) : 1.f) << ", ideal = " << m_ratio << ", nominal = " << getRatio() << ")" << std::endl;
588 // } else {
589 // std::cerr << " !transient (count = " << count << ", prev = " << m_prevTransientScore << ", diff = " << count - m_prevTransientScore << ")" << std::endl;
590 #endif
591 }
592
593 m_prevTransientScore = count;
594
595 return isTransient;
596 }
597
598 void
synthesiseBlock(size_t c,float * out,float * modulation,size_t lastStep)599 AudioTimeStretcher::synthesiseBlock(size_t c,
600 float *out,
601 float *modulation,
602 size_t lastStep)
603 {
604 bool unchanged = (lastStep == m_n1);
605
606 for (size_t i = 0; i <= m_wlen/2; ++i) {
607
608 float phase = princargf(atan2f(m_freq[c][i][1], m_freq[c][i][0]));
609 float adjustedPhase = phase;
610
611 // float binfreq = float(m_sampleRate * i) / m_wlen;
612
613 if (!unchanged) {
614
615 float mag = sqrtf(m_freq[c][i][0] * m_freq[c][i][0] +
616 m_freq[c][i][1] * m_freq[c][i][1]);
617
618 float omega = (2 * M_PI * m_n1 * i) / m_wlen;
619
620 float expectedPhase = m_prevPhase[c][i] + omega;
621
622 float phaseError = princargf(phase - expectedPhase);
623
624 float phaseIncrement = (omega + phaseError) / m_n1;
625
626 adjustedPhase = m_prevAdjustedPhase[c][i] +
627 lastStep * phaseIncrement;
628
629 float real = mag * cosf(adjustedPhase);
630 float imag = mag * sinf(adjustedPhase);
631 m_freq[c][i][0] = real;
632 m_freq[c][i][1] = imag;
633 }
634
635 m_prevPhase[c][i] = phase;
636 m_prevAdjustedPhase[c][i] = adjustedPhase;
637 }
638
639 fftwf_execute(m_iplan[c]); // m_freq -> m_time, inverse fft
640
641 for (size_t i = 0; i < m_wlen/2; ++i) {
642 float temp = m_time[c][i];
643 m_time[c][i] = m_time[c][i + m_wlen/2];
644 m_time[c][i + m_wlen/2] = temp;
645 }
646
647 for (size_t i = 0; i < m_wlen; ++i) {
648 m_time[c][i] = m_time[c][i] / m_wlen;
649 }
650
651 m_synthesisWindow->cut(m_time[c]);
652
653 for (size_t i = 0; i < m_wlen; ++i) {
654 out[i] += m_time[c][i];
655 }
656
657 if (modulation) {
658
659 float area = m_analysisWindow->getArea();
660
661 for (size_t i = 0; i < m_wlen; ++i) {
662 float val = m_synthesisWindow->getValue(i);
663 modulation[i] += val * area;
664 }
665 }
666 }
667
668
669
670 }
671
672