-
Notifications
You must be signed in to change notification settings - Fork 0
/
slant_edge.cpp
140 lines (116 loc) · 4.46 KB
/
slant_edge.cpp
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
// Measure the MTF from a slant-edge region.
// Author: Philip Salvaggio
#include "mats.h"
#include <iostream>
#include <gflags/gflags.h>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
DEFINE_string(image_file, "", "REQUIRED: Image filename.");
DEFINE_string(config_file, "", "Optional: SimulationConfig filename for a "
"plot overlay.");
DEFINE_int32(simulation_id, -1, "Simulation ID in config_file "
"(defaults to first).");
DEFINE_double(orientation, 0, "Orientation of the aperture.");
DEFINE_bool(whole_image, false, "Use the entire image.");
DEFINE_bool(quiet, false, "Do not plot");
DEFINE_bool(output_esf, false, "Output the ESF instead of the MTF.");
using namespace cv;
using namespace std;
int main(int argc, char** argv) {
google::ParseCommandLineFlags(&argc, &argv, true);
if (FLAGS_image_file == "") {
cerr << "Image filename is required." << endl;
return 1;
}
SlantEdgeMtf mtf_measurer;
Mat image = imread(mats::ResolvePath(FLAGS_image_file), 0);
if (!image.data) {
cerr << "Could not read image file." << endl;
return 1;
}
// Extract the slant edge ROI from the image.
Mat roi;
if (FLAGS_whole_image) {
roi = image;
} else {
auto bounds = move(GetRoi(image));
image(Range(bounds[1], bounds[1] + bounds[3]),
Range(bounds[0], bounds[0] + bounds[2])).copyTo(roi);
}
SlantEdgeMtf mtf_analyzer;
if (FLAGS_output_esf) {
double edge[3];
if (mtf_analyzer.DetectEdge(roi, edge)) {
int samples = mtf_analyzer.GetSamplesPerPixel(image, edge);
vector<double> esf, esf_stddevs;
mtf_analyzer.GenerateEsf(roi, edge, samples, &esf, &esf_stddevs);
mtf_analyzer.SmoothEsf(&esf);
for (size_t i = 0; i < esf.size(); i++) {
cout << i / double(samples) << "\t" << esf[i] << endl;
}
}
return 0;
}
// Perform the slant edge analysis on the image.
double orientation;
vector<double> mtf;
mtf_analyzer.Analyze(roi, &orientation, &mtf);
// Create the plot data and print out the MTF.
vector<pair<double, double>> mtf_data;
for (size_t i = 0; i < mtf.size(); i++) {
double freq = i / (2. * (mtf.size() - 1));
mtf_data.emplace_back(freq, mtf[i]);
cout << freq << "\t" << mtf[i] << endl;
}
if (FLAGS_quiet) return 0;
Gnuplot gp;
gp << "set xlabel \"Spatial Frequency [cyc/pixel]\"\n"
<< "set ylabel \"MTF\"\n";
if (FLAGS_config_file == "") {
gp << "unset key\n"
<< "plot " << gp.file1d(mtf_data) << " w l lw 3\n"
<< endl;
} else {
mats::SimulationConfig sim_config;
mats::DetectorParameters detector_params;
if (!mats::MatsInit(mats::ResolvePath(FLAGS_config_file),
&sim_config,
&detector_params)) {
return 1;
}
// Initialize the simulation parameters.
sim_config.set_array_size(512);
int sim_index = mats::LookupSimulationId(sim_config, FLAGS_simulation_id);
mats::Telescope telescope(sim_config, sim_index, detector_params);
telescope.detector()->set_rows(512);
telescope.detector()->set_cols(512);
// Set up the spectral resolution of the simulation.
vector<vector<double>> raw_weighting;
mats_io::TextFileReader::Parse(
sim_config.spectral_weighting_filename(),
&raw_weighting);
const vector<double>& wavelengths(raw_weighting[0]),
spectral_weighting(raw_weighting[1]);
double q = telescope.EffectiveQ(wavelengths, spectral_weighting);
cerr << "F/#: " << telescope.FNumber() << endl;
cerr << "Effective Q: " << q << endl;
Mat_<complex<double>> theoretical_otf;
telescope.EffectiveOtf(wavelengths, spectral_weighting, 0, 0,
&theoretical_otf);
Mat theoretical_2d_mtf = magnitude(theoretical_otf);
// Grab the radial profile of the OTF and convert to MTF.
std::vector<double> theoretical_mtf;
GetRadialProfile(FFTShift(theoretical_2d_mtf),
FLAGS_orientation * M_PI / 180,
&theoretical_mtf);
vector<pair<double, double>> t_mtf_data;
for (size_t i = 0; i < theoretical_mtf.size(); i++) {
double pix_freq = i / (2. * (theoretical_mtf.size() - 1)); // [cyc/pixel]
t_mtf_data.emplace_back(pix_freq, theoretical_mtf[i]);
}
gp << "plot " << gp.file1d(mtf_data) << " w l lw 3 t \"Measured\", "
<< gp.file1d(t_mtf_data) << " w l lw 3 t \"Theoretical\"\n"
<< endl;
}
return 0;
}