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 159 160 161
|
#ifndef NORMALIZE_PASSBAND_H
#define NORMALIZE_PASSBAND_H
#include "../algorithms/thresholdtools.h"
#include <aocommon/uvector.h>
#include <vector>
#ifdef __SSE__
#define USE_INTRINSICS
#endif
#ifdef USE_INTRINSICS
#include <xmmintrin.h>
#endif
class NormalizePassband
{
public:
static void NormalizeStepwise(TimeFrequencyData& data, size_t steps)
{
const size_t height = data.ImageHeight();
aocommon::UVector<num_t> stddev(steps);
for(size_t step=0; step!=steps; ++step)
{
const size_t startY = step*height/steps, endY = (step+1)*height/steps;
std::vector<num_t> dataVector((1+endY-startY) * data.ImageWidth() * data.ImageCount());
std::vector<num_t>::iterator vecIter = dataVector.begin();
const Mask2DCPtr mask = data.GetSingleMask();
for(size_t i=0; i!=data.ImageCount(); ++i)
{
const Image2D &image = *data.GetImage(i);
for(size_t y=startY; y!=endY; ++y)
{
const num_t* inputPtr = image.ValuePtr(0, y);
const bool* maskPtr = mask->ValuePtr(0, y);
for(size_t x=0; x!=image.Width(); ++x)
{
if(!*maskPtr && std::isfinite(*inputPtr))
{
*vecIter = *inputPtr;
++vecIter;
}
++inputPtr;
++maskPtr;
}
}
}
dataVector.resize(vecIter - dataVector.begin());
num_t mean;
ThresholdTools::WinsorizedMeanAndStdDev<num_t>(dataVector, mean, stddev[step]);
}
for(size_t i=0; i!=data.ImageCount(); ++i)
{
const Image2D &image = *data.GetImage(i);
Image2DPtr destImage = Image2D::CreateUnsetImagePtr(image.Width(), image.Height());
for(size_t step=0; step!=steps; ++step)
{
const size_t startY = step*height/steps, endY = (step+1)*height/steps;
float correctionFactor;
if(stddev[step] == 0.0)
correctionFactor = 0.0;
else
correctionFactor = 1.0 / stddev[step];
#ifdef USE_INTRINSICS
const __m128 corrFact4 = _mm_set_ps(correctionFactor, correctionFactor, correctionFactor, correctionFactor);
#endif
for(size_t y=startY; y!=endY; ++y)
{
const float *inputPtr = image.ValuePtr(0, y);
float *destPtr = destImage->ValuePtr(0, y);
#ifdef USE_INTRINSICS
for(size_t x=0;x<image.Width();x+=4)
{
_mm_store_ps(destPtr, _mm_mul_ps(corrFact4, _mm_load_ps(inputPtr)));
inputPtr += 4;
destPtr += 4;
}
#else
for(size_t x=0;x<image.Width();x++)
{
*destPtr = correctionFactor * *inputPtr;
inputPtr ++;
destPtr ++;
}
#endif
}
}
data.SetImage(i, std::move(destImage));
}
}
static void NormalizeSmooth(TimeFrequencyData& data)
{
TimeFrequencyData
real = data.Make(TimeFrequencyData::RealPart),
imag = data.Make(TimeFrequencyData::ImaginaryPart);
const size_t height = data.ImageHeight();
aocommon::UVector<double> rms(height);
std::vector<std::complex<num_t>> dataVector;
for(size_t y=0; y!=height; ++y)
{
dataVector.resize(data.ImageWidth() * data.ImageCount());
auto vecIter = dataVector.begin();
const Mask2DCPtr mask = data.GetSingleMask();
for(size_t i=0; i!=real.ImageCount(); ++i)
{
const Image2D& realImage = *real.GetImage(i);
const Image2D& imagImage = *imag.GetImage(i);
const num_t* realPtr = realImage.ValuePtr(0, y);
const num_t* imagPtr = imagImage.ValuePtr(0, y);
const bool* maskPtr = mask->ValuePtr(0, y);
for(size_t x=0; x!=realImage.Width(); ++x)
{
if(!*maskPtr && std::isfinite(*realPtr) && std::isfinite(*imagPtr))
{
*vecIter = std::complex<num_t>(*realPtr, *imagPtr);
++vecIter;
}
++realPtr;
++imagPtr;
++maskPtr;
}
}
dataVector.resize(vecIter - dataVector.begin());
rms[y] = ThresholdTools::WinsorizedRMS<num_t>(dataVector);
}
for(size_t i=0; i!=data.ImageCount(); ++i)
{
const Image2D &image = *data.GetImage(i);
Image2DPtr destImage = Image2D::CreateUnsetImagePtr(image.Width(), image.Height());
for(size_t y=0; y!=height; ++y)
{
num_t correctionFactor;
if(rms[y] == 0.0)
correctionFactor = 0.0;
else
correctionFactor = 1.0 / rms[y];
const num_t *inputPtr = image.ValuePtr(0, y);
num_t *destPtr = destImage->ValuePtr(0, y);
for(size_t x=0;x<image.Width();x++)
{
*destPtr = correctionFactor * *inputPtr;
++inputPtr;
++destPtr;
}
}
data.SetImage(i, std::move(destImage));
}
}
};
#endif
|