File: resampling.h

package info (click to toggle)
aoflagger 3.4.0-4
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 8,960 kB
  • sloc: cpp: 83,076; python: 10,187; sh: 260; makefile: 178
file content (158 lines) | stat: -rw-r--r-- 5,926 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
#ifndef AOFLAGGER_ALGORITHMS_RESAMPLING_H
#define AOFLAGGER_ALGORITHMS_RESAMPLING_H

#include <algorithm>
#include <utility>
#include <vector>

#include "../structures/timefrequencydata.h"
#include "../structures/timefrequencymetadata.h"

#include "thresholdtools.h"

namespace algorithms {

void downsample_masked(TimeFrequencyData& tfData,
                       TimeFrequencyMetaData* metaData, size_t horizontalFactor,
                       size_t verticalFactor) {
  // Decrease in horizontal direction
  size_t polCount = tfData.PolarizationCount();
  for (size_t i = 0; i < polCount; ++i) {
    TimeFrequencyData polData(tfData.MakeFromPolarizationIndex(i));
    const Mask2DCPtr mask = polData.GetSingleMask();
    for (unsigned j = 0; j < polData.ImageCount(); ++j) {
      const Image2DCPtr image = polData.GetImage(j);
      polData.SetImage(j, ThresholdTools::ShrinkHorizontally(
                              horizontalFactor, image.get(), mask.get()));
    }
    tfData.SetPolarizationData(i, std::move(polData));
  }
  size_t maskCount = tfData.MaskCount();
  for (size_t i = 0; i < maskCount; ++i) {
    Mask2DCPtr mask = tfData.GetMask(i);
    Mask2DPtr newMask(
        new Mask2D(mask->ShrinkHorizontallyForAveraging(horizontalFactor)));
    tfData.SetMask(i, std::move(newMask));
  }

  // Decrease in vertical direction
  for (size_t i = 0; i < polCount; ++i) {
    TimeFrequencyData polData(tfData.MakeFromPolarizationIndex(i));
    const Mask2DCPtr mask = polData.GetSingleMask();
    for (unsigned j = 0; j < polData.ImageCount(); ++j) {
      const Image2DCPtr image = polData.GetImage(j);
      polData.SetImage(j, ThresholdTools::ShrinkVertically(
                              verticalFactor, image.get(), mask.get()));
    }
    tfData.SetPolarizationData(i, std::move(polData));
  }
  for (size_t i = 0; i < maskCount; ++i) {
    Mask2DCPtr mask = tfData.GetMask(i);
    Mask2DPtr newMask(
        new Mask2D(mask->ShrinkVerticallyForAveraging(verticalFactor)));
    tfData.SetMask(i, std::move(newMask));
  }

  if (metaData) {
    if (metaData->HasBand() && verticalFactor != 1) {
      BandInfo newBand = metaData->Band();
      size_t newNChannels = tfData.ImageHeight();
      newBand.channels.resize(newNChannels);
      for (size_t i = 0; i != newNChannels; ++i) {
        const size_t startChannel = i * verticalFactor;
        const size_t endChannel = std::min((i + 1) * verticalFactor,
                                           metaData->Band().channels.size());
        double f = 0;
        for (size_t j = startChannel; j != endChannel; ++j) {
          f += metaData->Band().channels[j].frequencyHz;
        }
        newBand.channels[i].frequencyHz = f / (endChannel - startChannel);
      }
      metaData->SetBand(newBand);
    }
    if (metaData->HasObservationTimes() && horizontalFactor != 1) {
      size_t newNTimes = tfData.ImageWidth();
      std::vector<double> times(newNTimes);
      for (size_t i = 0; i != newNTimes; ++i) {
        const size_t startTime = i * horizontalFactor;
        const size_t endTime = std::min((i + 1) * horizontalFactor,
                                        metaData->ObservationTimes().size());
        double t = 0;
        for (size_t j = startTime; j != endTime; ++j) {
          t += metaData->ObservationTimes()[j];
        }
        times[i] = t / (endTime - startTime);
      }
      metaData->SetObservationTimes(times);
    }
  }
}

void upsample_image(TimeFrequencyData timeFrequencyData,
                    TimeFrequencyData& destination, size_t horizontalFactor,
                    size_t verticalFactor) {
  const size_t imageCount = timeFrequencyData.ImageCount(),
               newWidth = destination.ImageWidth(),
               newHeight = destination.ImageHeight();
  if (destination.ImageCount() != imageCount)
    throw std::runtime_error(
        "Error in upsample() call: source and image have different number of "
        "images");

  if (horizontalFactor > 1) {
    for (size_t i = 0; i < imageCount; ++i) {
      Image2DPtr newImage(
          new Image2D(timeFrequencyData.GetImage(i)->EnlargeHorizontally(
              horizontalFactor, newWidth)));
      timeFrequencyData.SetImage(i, newImage);
    }
  }

  for (size_t i = 0; i < imageCount; ++i) {
    if (verticalFactor > 1) {
      Image2DPtr newImage(
          new Image2D(timeFrequencyData.GetImage(i)->EnlargeVertically(
              verticalFactor, newHeight)));
      destination.SetImage(i, newImage);
    } else {
      destination.SetImage(i, timeFrequencyData.GetImage(i));
    }
  }
}

void upsample_mask(TimeFrequencyData timeFrequencyData,
                   TimeFrequencyData& destination, size_t horizontalFactor,
                   size_t verticalFactor) {
  const size_t maskCount = timeFrequencyData.MaskCount(),
               newWidth = destination.ImageWidth(),
               newHeight = destination.ImageHeight(),
               oldHeight = timeFrequencyData.ImageHeight();
  if (destination.MaskCount() != maskCount)
    throw std::runtime_error(
        "Error in upsample() call: source and image have different number of "
        "masks");

  if (horizontalFactor > 1) {
    for (size_t i = 0; i != maskCount; ++i) {
      Mask2DPtr newMask = Mask2D::CreateUnsetMaskPtr(newWidth, oldHeight);
      newMask->EnlargeHorizontallyAndSet(*timeFrequencyData.GetMask(i),
                                         horizontalFactor);
      timeFrequencyData.SetMask(i, newMask);
    }
  }

  for (size_t i = 0; i != maskCount; ++i) {
    if (verticalFactor > 1) {
      Mask2DPtr newMask = Mask2D::CreateUnsetMaskPtr(newWidth, newHeight);
      newMask->EnlargeVerticallyAndSet(*timeFrequencyData.GetMask(i),
                                       verticalFactor);
      destination.SetMask(i, newMask);
    } else {
      destination.SetMask(i, timeFrequencyData.GetMask(i));
    }
  }
}

}  // namespace algorithms

#endif