← Back to branch summary

~ubuntu-branches/ubuntu/wily/openms/wily

« back to all changes in this revision

Viewing changes to source/ANALYSIS/MAPMATCHING/PoseClusteringShiftSuperimposer.C

  • Committer: Package Import Robot
  • Author(s): Filippo Rusconi
  • Date: 2013-12-20 11:30:16 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: package-import@ubuntu.com-20131220113016-wre5g9bteeheq6he
Tags: 1.11.1-3
* remove version number from libbost development package names;
* ensure that AUTHORS is correctly shipped in all packages.

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/ANALYSIS/MAPMATCHING/PoseClusteringShiftSuperimposer.C
1
 
// -*- mode: C++; tab-width: 2; -*-
2
 
// vi: set ts=2:
3
 
//
4
 
// --------------------------------------------------------------------------
5
 
//                   OpenMS Mass Spectrometry Framework
6
 
// --------------------------------------------------------------------------
7
 
//  Copyright (C) 2003-2011 -- Oliver Kohlbacher, Knut Reinert
8
 
//
9
 
//  This library is free software; you can redistribute it and/or
10
 
//  modify it under the terms of the GNU Lesser General Public
11
 
//  License as published by the Free Software Foundation; either
12
 
//  version 2.1 of the License, or (at your option) any later version.
13
 
//
14
 
//  This library is distributed in the hope that it will be useful,
15
 
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
16
 
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17
 
//  Lesser General Public License for more details.
18
 
//
19
 
//  You should have received a copy of the GNU Lesser General Public
20
 
//  License along with this library; if not, write to the Free Software
21
 
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
1
// --------------------------------------------------------------------------
 
2
//                   OpenMS -- Open-Source Mass Spectrometry
 
3
// --------------------------------------------------------------------------
 
4
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
 
5
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
 
6
//
 
7
// This software is released under a three-clause BSD license:
 
8
//  * Redistributions of source code must retain the above copyright
 
9
//    notice, this list of conditions and the following disclaimer.
 
10
//  * Redistributions in binary form must reproduce the above copyright
 
11
//    notice, this list of conditions and the following disclaimer in the
 
12
//    documentation and/or other materials provided with the distribution.
 
13
//  * Neither the name of any author or any participating institution
 
14
//    may be used to endorse or promote products derived from this software
 
15
//    without specific prior written permission.
 
16
// For a full list of authors, refer to the file AUTHORS.
 
17
// --------------------------------------------------------------------------
 
18
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 
19
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
20
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
21
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
 
22
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 
23
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 
24
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 
25
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 
26
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 
27
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 
28
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
22
29
//
23
30
// --------------------------------------------------------------------------
24
31
// $Maintainer: Clemens Groepl $
42
49
#else
43
50
#define V_(bla)
44
51
#endif
45
 
#define VV_(bla) V_(""#bla": " << bla)
 
52
#define VV_(bla) V_("" # bla ": " << bla)
46
53
 
47
54
namespace OpenMS
48
55
{
53
60
    setName(getProductName());
54
61
 
55
62
    defaults_.setValue("mz_pair_max_distance", 0.5, "Maximum of m/z deviation of corresponding elements in different maps.  "
56
 
      "This condition applies to the pairs considered in hashing.");
 
63
                                                    "This condition applies to the pairs considered in hashing.");
57
64
    defaults_.setMinFloat("mz_pair_max_distance", 0.);
58
65
 
59
66
    defaults_.setValue("num_used_points", 2000, "Maximum number of elements considered in each map "
60
 
      "(selected by intensity).  Use this to reduce the running time "
61
 
      "and to disregard weak signals during alignment.  For using all points, set this to -1.");
 
67
                                                "(selected by intensity).  Use this to reduce the running time "
 
68
                                                "and to disregard weak signals during alignment.  For using all points, set this to -1.");
62
69
    defaults_.setMinInt("num_used_points", -1);
63
70
 
64
71
    defaults_.setValue("shift_bucket_size", 3.0, "The shift of the retention time "
65
 
      "interval is being hashed into buckets of this size during pose "
66
 
      "clustering.  A good choice for this would be about "
67
 
      "the time between consecutive MS scans.");
 
72
                                                 "interval is being hashed into buckets of this size during pose "
 
73
                                                 "clustering.  A good choice for this would be about "
 
74
                                                 "the time between consecutive MS scans.");
68
75
    defaults_.setMinFloat("shift_bucket_size", 0.);
69
76
 
70
77
    defaults_.setValue("max_shift", 1000.0, "Maximal shift which is considered during histogramming.  "
71
 
      "This applies for both directions.", StringList::create("advanced"));
 
78
                                            "This applies for both directions.", StringList::create("advanced"));
72
79
    defaults_.setMinFloat("max_shift", 0.);
73
80
 
74
81
    defaults_.setValue("dump_buckets", "", "[DEBUG] If non-empty, base filename where hash table buckets will be dumped to.  "
75
 
      "A serial number for each invocation will be appended automatically.", StringList::create("advanced"));
 
82
                                           "A serial number for each invocation will be appended automatically.", StringList::create("advanced"));
76
83
 
77
84
    defaults_.setValue("dump_pairs", "", "[DEBUG] If non-empty, base filename where the individual hashed pairs will be dumped to (large!).  "
78
 
      "A serial number for each invocation will be appended automatically.", StringList::create("advanced"));
 
85
                                         "A serial number for each invocation will be appended automatically.", StringList::create("advanced"));
79
86
 
80
87
    defaultsToParam_();
81
88
    return;
82
89
  }
83
90
 
84
 
  void
85
 
  PoseClusteringShiftSuperimposer::run(const std::vector<ConsensusMap>& maps, std::vector<TransformationDescription>& transformations)
 
91
  void PoseClusteringShiftSuperimposer::run(const ConsensusMap & map_model, const ConsensusMap & map_scene, TransformationDescription & transformation)
86
92
  {
87
 
 
88
 
    if ( maps.size() != 2 )
89
 
    {
90
 
      throw Exception::IllegalArgument(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Exactly two input maps are required");
91
 
    }
92
 
 
93
93
    typedef ConstRefVector<ConsensusMap> PeakPointerArray_;
94
94
    typedef Math::LinearInterpolation<DoubleReal, DoubleReal> LinearInterpolationType_;
95
95
 
119
119
    // Optionally, we will write dumps of the hash table buckets.
120
120
    bool do_dump_buckets = false;
121
121
    String dump_buckets_basename;
122
 
    if ( param_.getValue("dump_buckets") != "" )
 
122
    if (param_.getValue("dump_buckets") != "")
123
123
    {
124
124
      do_dump_buckets = true;
125
125
      dump_buckets_basename = param_.getValue("dump_buckets");
129
129
    // Even more optionally, we will write dumps of the hashed pairs.
130
130
    bool do_dump_pairs = false;
131
131
    String dump_pairs_basename;
132
 
    if ( param_.getValue("dump_pairs") != "" )
 
132
    if (param_.getValue("dump_pairs") != "")
133
133
    {
134
134
      do_dump_pairs = true;
135
135
      dump_pairs_basename = param_.getValue("dump_pairs");
139
139
    //**************************************************************************
140
140
    // Select the most abundant data points only.  After that, disallow modifications
141
141
    // (we tend to have annoying issues with const_iterator versus iterator).
142
 
    PeakPointerArray_ model_map_ini(maps[0].begin(), maps[0].end());
 
142
    PeakPointerArray_ model_map_ini(map_model.begin(), map_model.end());
143
143
    const PeakPointerArray_ & model_map(model_map_ini);
144
 
    PeakPointerArray_ scene_map_ini(maps[1].begin(), maps[1].end());
 
144
    PeakPointerArray_ scene_map_ini(map_scene.begin(), map_scene.end());
145
145
    const PeakPointerArray_ & scene_map(scene_map_ini);
146
146
    {
147
147
      // truncate the data as necessary
148
148
      // casting to SignedSize is done on PURPOSE here! (num_used_points will be maximal if -1 is used)
149
149
      const Size num_used_points = (SignedSize) param_.getValue("num_used_points");
150
 
      if ( model_map_ini.size() > num_used_points )
 
150
      if (model_map_ini.size() > num_used_points)
151
151
      {
152
152
        model_map_ini.sortByIntensity(true);
153
153
        model_map_ini.resize(num_used_points);
154
154
      }
155
155
      model_map_ini.sortByComparator(Peak2D::MZLess());
156
156
      setProgress(++actual_progress);
157
 
      if ( scene_map_ini.size() > num_used_points )
 
157
      if (scene_map_ini.size() > num_used_points)
158
158
      {
159
159
        scene_map_ini.sortByIntensity(true);
160
160
        scene_map_ini.resize(num_used_points);
171
171
    // get RT ranges (NOTE: we trust that min and max have been updated in the
172
172
    // ConsensusMap::convert() method !)
173
173
 
174
 
    const DoubleReal model_low = maps[0].getMin()[ConsensusFeature::RT];
175
 
    const DoubleReal scene_low = maps[1].getMin()[ConsensusFeature::RT];
176
 
    const DoubleReal model_high = maps[0].getMax()[ConsensusFeature::RT];
177
 
    const DoubleReal scene_high = maps[1].getMax()[ConsensusFeature::RT];
 
174
    const DoubleReal model_low = map_model.getMin()[ConsensusFeature::RT];
 
175
    const DoubleReal scene_low = map_scene.getMin()[ConsensusFeature::RT];
 
176
    const DoubleReal model_high = map_model.getMax()[ConsensusFeature::RT];
 
177
    const DoubleReal scene_high = map_scene.getMax()[ConsensusFeature::RT];
178
178
 
179
179
    // OLD STUFF
180
180
    //    const DoubleReal rt_low = (maps[0].getMin()[ConsensusFeature::RT] + maps[1].getMin()[ConsensusFeature::RT]) / 2.;
188
188
      // actually the largest possible shift can be much smaller, depending on the data
189
189
      do
190
190
      {
191
 
        if ( max_shift < 0 )
 
191
        if (max_shift < 0)
192
192
          max_shift = -max_shift;
193
193
        //     ...ml@@@mh........    ,    ........ml@@@mh...
194
194
        //     ........sl@@@sh...    ,    ...sl@@@sh........
195
195
        DoubleReal diff;
196
196
        diff = model_high - scene_low;
197
 
        if ( diff < 0 )
 
197
        if (diff < 0)
198
198
          diff = -diff;
199
 
        if ( max_shift > diff )
 
199
        if (max_shift > diff)
200
200
          max_shift = diff;
201
201
        diff = model_low - scene_high;
202
 
        if ( diff < 0 )
 
202
        if (diff < 0)
203
203
          diff = -diff;
204
 
        if ( max_shift > diff )
 
204
        if (max_shift > diff)
205
205
          max_shift = diff;
206
206
      }
207
 
      while ( 0 );
 
207
      while (0);
208
208
 
209
209
      const Int shift_buckets_num_half = 4 + (Int) ceil((max_shift) / shift_bucket_size);
210
210
      const Int shift_buckets_num = 1 + 2 * shift_buckets_num_half;
221
221
    do
222
222
    {
223
223
      DoubleReal total_int_model_map = 0;
224
 
      for ( Size i = 0; i < model_map.size(); ++i )
 
224
      for (Size i = 0; i < model_map.size(); ++i)
225
225
      {
226
226
        total_int_model_map += model_map[i].getIntensity();
227
227
      }
228
228
      setProgress(++actual_progress);
229
229
      DoubleReal total_int_scene_map = 0;
230
 
      for ( Size i = 0; i < scene_map.size(); ++i )
 
230
      for (Size i = 0; i < scene_map.size(); ++i)
231
231
      {
232
232
        total_int_scene_map += scene_map[i].getIntensity();
233
233
      }
235
235
      // ... and finally ...
236
236
      total_intensity_ratio = total_int_model_map / total_int_scene_map;
237
237
    }
238
 
    while ( 0 ); // (the extra syntax helps with code folding in eclipse!)
 
238
    while (0);   // (the extra syntax helps with code folding in eclipse!)
239
239
    setProgress((actual_progress = 20));
240
240
 
241
241
    /// The serial number is incremented for each invocation of this, to avoid overwriting of hash table dumps.
266
266
    {
267
267
      String dump_pairs_filename;
268
268
      std::ofstream dump_pairs_file;
269
 
      if ( do_dump_pairs )
 
269
      if (do_dump_pairs)
270
270
      {
271
271
        dump_pairs_filename = dump_pairs_basename + String(dump_buckets_serial);
272
272
        dump_pairs_file.open(dump_pairs_filename.c_str());
275
275
      setProgress(++actual_progress);
276
276
 
277
277
      // first point in model map
278
 
      for ( Size i = 0, i_low = 0, i_high = 0, k_low = 0, k_high = 0; i < model_map_size - 1; ++i )
 
278
      for (Size i = 0, i_low = 0, i_high = 0, k_low = 0, k_high = 0; i < model_map_size - 1; ++i)
279
279
      {
280
280
        setProgress(actual_progress + Real(i) / model_map_size * 10.f);
281
281
 
282
282
        // Adjust window around i in model map
283
 
        while ( i_low < model_map_size && model_map[i_low].getMZ() < model_map[i].getMZ() - mz_pair_max_distance )
 
283
        while (i_low < model_map_size && model_map[i_low].getMZ() < model_map[i].getMZ() - mz_pair_max_distance)
284
284
          ++i_low;
285
 
        while ( i_high < model_map_size && model_map[i_high].getMZ() <= model_map[i].getMZ() + mz_pair_max_distance )
 
285
        while (i_high < model_map_size && model_map[i_high].getMZ() <= model_map[i].getMZ() + mz_pair_max_distance)
286
286
          ++i_high;
287
287
        DoubleReal i_winlength_factor = 1. / (i_high - i_low);
288
288
        i_winlength_factor -= winlength_factor_baseline;
289
 
        if ( i_winlength_factor <= 0 )
 
289
        if (i_winlength_factor <= 0)
290
290
          continue;
291
291
 
292
292
        // Adjust window around k in scene map
293
 
        while ( k_low < scene_map_size && scene_map[k_low].getMZ() < model_map[i].getMZ() - mz_pair_max_distance )
 
293
        while (k_low < scene_map_size && scene_map[k_low].getMZ() < model_map[i].getMZ() - mz_pair_max_distance)
294
294
          ++k_low;
295
 
        while ( k_high < scene_map_size && scene_map[k_high].getMZ() <= model_map[i].getMZ() + mz_pair_max_distance )
 
295
        while (k_high < scene_map_size && scene_map[k_high].getMZ() <= model_map[i].getMZ() + mz_pair_max_distance)
296
296
          ++k_high;
297
297
 
298
298
        // first point in scene map
299
 
        for ( Size k = k_low; k < k_high; ++k )
 
299
        for (Size k = k_low; k < k_high; ++k)
300
300
        {
301
301
          DoubleReal k_winlength_factor = 1. / (k_high - k_low);
302
302
          k_winlength_factor -= winlength_factor_baseline;
303
 
          if ( k_winlength_factor <= 0 )
 
303
          if (k_winlength_factor <= 0)
304
304
            continue;
305
305
 
306
306
          // compute similarity of intensities i k
321
321
          // hash the images of scaling, rt_low and rt_high into their respective hash tables
322
322
          shift_hash_.addValue(shift, similarity_ik);
323
323
 
324
 
          if ( do_dump_pairs )
 
324
          if (do_dump_pairs)
325
325
          {
326
326
            dump_pairs_file << i << ' ' << model_map[i].getRT() << ' ' << model_map[i].getMZ() << ' ' << k << ' ' << scene_map[k].getRT() << ' '
327
 
                << scene_map[k].getMZ() << ' ' << similarity_ik << ' ' << std::endl;
 
327
                            << scene_map[k].getMZ() << ' ' << similarity_ik << ' ' << std::endl;
328
328
          }
329
329
 
330
330
        } // k
331
331
      } // i
332
332
    }
333
 
    while ( 0 ); // end of hashing (the extra syntax helps with code folding in eclipse!)
 
333
    while (0);   // end of hashing (the extra syntax helps with code folding in eclipse!)
334
334
 
335
335
    setProgress((actual_progress = 30));
336
336
 
337
337
    ///////////////////////////////////////////////////////////////////
338
338
    // work on shift_hash_
339
 
//    DoubleReal shift_low;
340
 
//    DoubleReal shift_centroid;
341
 
//    DoubleReal shift_high;
 
339
    //   DoubleReal shift_low;
 
340
    //   DoubleReal shift_centroid;
 
341
    //   DoubleReal shift_high;
342
342
 
343
 
// OLD STUFF
 
343
    // OLD STUFF
344
344
    DoubleReal shift_low;
345
345
    DoubleReal shift_centroid;
346
346
    DoubleReal shift_high;
352
352
      // optionally, dump before filtering
353
353
      String dump_buckets_filename;
354
354
      std::ofstream dump_buckets_file;
355
 
      if ( do_dump_buckets )
 
355
      if (do_dump_buckets)
356
356
      {
357
357
        dump_buckets_filename = dump_buckets_basename + "_" + String(dump_buckets_serial);
358
358
        dump_buckets_file.open(dump_buckets_filename.c_str());
360
360
 
361
361
        dump_buckets_file << "# shift hash table buckets dump ( scale, height ) : " << dump_buckets_filename << std::endl;
362
362
        dump_buckets_file << "# unfiltered hash data\n";
363
 
        for ( Size index = 0; index < shift_hash_.getData().size(); ++index )
 
363
        for (Size index = 0; index < shift_hash_.getData().size(); ++index)
364
364
        {
365
365
          const DoubleReal image = shift_hash_.index2key(index);
366
366
          const DoubleReal height = shift_hash_.getData()[index];
385
385
      shift_hash_.getData().swap(buffer);
386
386
 
387
387
      // optionally, dump after filtering
388
 
      if ( do_dump_buckets )
 
388
      if (do_dump_buckets)
389
389
      {
390
390
        dump_buckets_file << "# tophat filtered hash data\n";
391
 
        for ( Size index = 0; index < shift_hash_.getData().size(); ++index )
 
391
        for (Size index = 0; index < shift_hash_.getData().size(); ++index)
392
392
        {
393
393
          const DoubleReal image = shift_hash_.index2key(index);
394
394
          const DoubleReal height = shift_hash_.getData()[index];
409
409
          std::sort(buffer.begin(), buffer.end(), std::greater<DoubleReal>());
410
410
          DoubleReal freq_intercept = shift_hash_.getData().front();
411
411
          DoubleReal freq_slope = (shift_hash_.getData().back() - shift_hash_.getData().front()) / DoubleReal(buffer.size())
412
 
              / scaling_histogram_crossing_slope;
413
 
          if ( !freq_slope || !buffer.size() )
 
412
                                  / scaling_histogram_crossing_slope;
 
413
          if (!freq_slope || !buffer.size())
414
414
          {
415
415
            // in fact these conditions are actually impossible, but let's be really sure ;-)
416
416
            freq_cutoff_low = 0;
418
418
          else
419
419
          {
420
420
            Size index = 1; // not 0 (!)
421
 
            while ( buffer[index] >= freq_intercept + freq_slope * DoubleReal(index) )
 
421
            while (buffer[index] >= freq_intercept + freq_slope * DoubleReal(index))
422
422
            {
423
423
              ++index;
424
424
            }
426
426
          }
427
427
        }
428
428
      }
429
 
      while ( 0 );
 
429
      while (0);
430
430
      setProgress(++actual_progress);
431
431
 
432
432
      // apply freq_cutoff, setting smaller values to zero
433
 
      for ( Size index = 0; index < shift_hash_.getData().size(); ++index )
 
433
      for (Size index = 0; index < shift_hash_.getData().size(); ++index)
434
434
      {
435
 
        if ( shift_hash_.getData()[index] < freq_cutoff_low )
 
435
        if (shift_hash_.getData()[index] < freq_cutoff_low)
436
436
        {
437
437
          shift_hash_.getData()[index] = 0;
438
438
        }
440
440
      setProgress(++actual_progress);
441
441
 
442
442
      // optionally, dump after noise filtering using freq_cutoff
443
 
      if ( do_dump_buckets )
 
443
      if (do_dump_buckets)
444
444
      {
445
445
        dump_buckets_file << "# after freq_cutoff, which is: " << freq_cutoff_low << '\n';
446
 
        for ( Size index = 0; index < shift_hash_.getData().size(); ++index )
 
446
        for (Size index = 0; index < shift_hash_.getData().size(); ++index)
447
447
        {
448
448
          const DoubleReal image = shift_hash_.index2key(index);
449
449
          const DoubleReal height = shift_hash_.getData()[index];
462
462
        Size data_range_end = data_size;
463
463
        DoubleReal mean;
464
464
        DoubleReal stdev;
465
 
        for ( UInt loop = 0; loop < loops_mean_stdev_cutoff; ++loop ) // MAGIC ALERT: number of loops
 
465
        for (UInt loop = 0; loop < loops_mean_stdev_cutoff; ++loop)   // MAGIC ALERT: number of loops
466
466
        {
467
467
          statistics.update(data_begin + data_range_begin, data_begin + data_range_end);
468
468
          mean = statistics.mean() + data_range_begin;
474
474
          shift_low = (outside_mean - outside_stdev);
475
475
          shift_centroid = (outside_mean);
476
476
          shift_high = (outside_mean + outside_stdev);
477
 
          if ( do_dump_buckets )
 
477
          if (do_dump_buckets)
478
478
          {
479
 
            dump_buckets_file << "# loop: " << loop << "  mean: " << outside_mean << "  stdev: " << outside_stdev << "  (mean-stdev): " << outside_mean
480
 
                - outside_stdev << "  (mean+stdev): " << outside_mean + outside_stdev << "  data_range_begin: " << data_range_begin << "  data_range_end: "
481
 
                << data_range_end << std::endl;
 
479
            dump_buckets_file << "# loop: " << loop << "  mean: " << outside_mean << "  stdev: " << outside_stdev << "  (mean-stdev): "
 
480
                              << outside_mean - outside_stdev << "  (mean+stdev): " << outside_mean + outside_stdev
 
481
                              << "  data_range_begin: " << data_range_begin << "  data_range_end: "
 
482
                              << data_range_end << std::endl;
482
483
          }
483
484
        }
484
485
        setProgress(++actual_progress);
485
486
      }
486
 
      if ( do_dump_buckets )
 
487
      if (do_dump_buckets)
487
488
      {
488
489
        dump_buckets_file << "# EOF" << std::endl;
489
490
        dump_buckets_file.close();
491
492
      setProgress(80);
492
493
 
493
494
    }
494
 
    while ( 0 );
 
495
    while (0);
495
496
 
496
497
    //************************************************************************************
497
498
    // Estimate transform
502
503
    intercept = shift_centroid;
503
504
#else // ooh, use maximum bins instead (Note: this is a fossil which would disregard most of the above computations!  The code is left here for developers/debugging only.)
504
505
    const Size rt_low_max_index = std::distance(shift_hash_.getData().begin(),
505
 
        std::max_element(shift_hash_.getData().begin(), shift_hash_.getData().end()));
 
506
                                                std::max_element(shift_hash_.getData().begin(), shift_hash_.getData().end()));
506
507
    intercept = shift_hash_.index2key(rt_low_max_index);
507
508
#endif
508
509
 
512
513
 
513
514
    // set trafo
514
515
    {
515
 
      transformations.clear();
516
 
 
517
 
                        Param params;
 
516
      Param params;
518
517
      params.setValue("slope", 1.0);
519
518
      params.setValue("intercept", intercept);
520
519
 
521
 
                        TransformationDescription trafo;
522
 
                        trafo.fitModel("linear", params);
523
 
                        transformations.push_back(trafo);
524
 
     }
 
520
      TransformationDescription trafo;
 
521
      trafo.fitModel("linear", params);
 
522
      transformation = trafo;
 
523
    }
525
524
 
526
525
    setProgress(++actual_progress);
527
 
 
528
526
    endProgress();
529
527
 
530
528
    return;