Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Complex value support #47

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
72 changes: 56 additions & 16 deletions cmd/mrpeek.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,14 @@ using namespace VT;

vector<std::string> colourmap_choices_std;
vector<const char*> colourmap_choices_cstr;
const char* complex_choices [] = {
"real",
"imag",
"abs",
"arg",
"colour",
nullptr
};

enum ArrowMode { ARROW_SLICEVOL, ARROW_COLOUR, ARROW_CROSSHAIR, N_ARROW_MODES };

Expand All @@ -40,6 +48,7 @@ void usage ()
++entry;
} while (entry->name && !entry->special && !entry->is_colour);
colourmap_choices_cstr.reserve (colourmap_choices_std.size() + 1);

for (const auto& s : colourmap_choices_std)
colourmap_choices_cstr.push_back (s.c_str());
colourmap_choices_cstr.push_back (nullptr);
Expand Down Expand Up @@ -103,6 +112,11 @@ void usage ()
". Default is " + colourmap_choices_std[0] + ".")
+ Argument ("name").type_choice (colourmap_choices_cstr.data())

+ Option ("complex",
"for complex numbers, select which aspect to show. Choices are: " + join (complex_choices, ",") +
". Default is " + complex_choices[0] + ".")
+ Argument ("type").type_choice (complex_choices)

+ Option ("focus",
"set focus (crosshairs) at specified position, as a comma-separated "
"list of values. Use empty entries to leave as default (e.g. '-focus ,,100' "
Expand Down Expand Up @@ -130,7 +144,8 @@ void usage ()


using value_type = float;
using ImageType = Image<value_type>;
using cpx_value_type = cfloat;
using ImageType = Image<cpx_value_type>;
using Reslicer = Adapter::Reslice<Interp::Nearest, ImageType>;
using LinearReslicer = Adapter::Reslice<Interp::Linear, ImageType>;
using CubicReslicer = Adapter::Reslice<Interp::Cubic, ImageType>;
Expand All @@ -142,7 +157,8 @@ using CubicReslicer = Adapter::Reslice<Interp::Cubic, ImageType>;
// These will need to be moved into a struct/class eventually...
int levels = 32;
int x_axis, y_axis, slice_axis = 2, plot_axis = slice_axis, vol_axis = -1;
value_type pmin = DEFAULT_PMIN, pmax = DEFAULT_PMAX, zoom = 1.0;
int display_type = 0;
default_type pmin = DEFAULT_PMIN, pmax = DEFAULT_PMAX, zoom = 1.0;
bool crosshair = true, colorbar = true, orthoview = true, interactive = true;
bool do_plot = false, show_image = true, interpolate = false, show_text = true;
vector<int> focus (3, 0); // relative to original image grid
Expand All @@ -151,6 +167,18 @@ Sixel::ColourMaps colourmaps;
Sixel::ColourMaps plot_cmaps;



inline value_type get_scalar (const cpx_value_type val) {
switch (display_type) {
case 0: return std::real (val);
case 1: return std::imag (val);
case 2: return std::abs (val);
case 3: return std::arg (val);
case 4: return std::abs (val);
}
return 0.0;
}

inline std::string move_down (int n) {
if (interactive)
return move_cursor (Down, n);
Expand All @@ -164,7 +192,7 @@ inline std::string move_down (int n) {
// calculate percentile of a list of numbers
// implementation based on `mrthreshold` - can be merged with Math::median in due course
template <class Container>
value_type percentile (Container& data, default_type percentile)
default_type percentile (Container& data, default_type percentile)
{
// ignore nan
auto isnotfinite = [](typename Container::value_type val) { return !std::isfinite(val); };
Expand Down Expand Up @@ -236,7 +264,15 @@ inline std::string show_focus (ImageType& image)
}
out += "] ";

out += "| value: " + str(image.value());
out += "| value [";
switch (display_type) {
case 0: out += "real"; break;
case 1: out += "imag"; break;
case 2: out += "abs"; break;
case 3: out += "arg"; break;
case 4: out += "colour"; break;
}
out+= "]: " + str(image.value());

return out;
}
Expand All @@ -249,7 +285,7 @@ inline Reslicer get_regridder (ImageType& image, int with_slice_axis)
Header header_target (image);
default_type original_extent;
for (int d = 0; d < 3; ++d) {
const float new_voxel_size = (d == with_slice_axis) ? image.spacing(d) : 1.0f/zoom;
const default_type new_voxel_size = (d == with_slice_axis) ? image.spacing(d) : 1.0f/zoom;

original_extent = image.size(d) * image.spacing(d);

Expand All @@ -276,16 +312,16 @@ void autoscale (ImageType& image, Sixel::CMap& cmap)
const int y_dim = image_regrid.size(y_axis);
image_regrid.index(slice_axis) = focus[slice_axis];

vector<value_type> currentslice (x_dim*y_dim);
vector<default_type> currentslice (x_dim*y_dim);
size_t k = 0;
std::cerr.flush();
for (auto l = Loop (vector<size_t>({ size_t(x_axis), size_t(y_axis) }))(image_regrid); l; ++l) {
//std::cerr << "[" << image_regrid.index(0) << " " << image_regrid.index(1) << " " << image_regrid.index(2) << "] ";
currentslice[k++] = image_regrid.value();
currentslice[k++] = get_scalar (image_regrid.value());
}

value_type vmin = percentile(currentslice, pmin);
value_type vmax = percentile(currentslice, pmax);
default_type vmin = percentile(currentslice, pmin);
default_type vmax = percentile(currentslice, pmax);
cmap.set_scaling_min_max (vmin, vmax);
INFO("reset intensity range to " + str(vmin) + " - " +str(vmax));
}
Expand Down Expand Up @@ -333,7 +369,7 @@ void render_slice (InterpType& regrid, const Sixel::ViewPort& view, const Sixel:
regrid.index(y_axis) = y_dim-1-y;
for (int x = 0; x < x_dim; ++x) {
regrid.index(x_axis) = x_dim-1-x;
view(x,y) = cmap (regrid.value());
view(x,y) = cmap (get_scalar (regrid.value()));
}
}
}
Expand Down Expand Up @@ -381,16 +417,16 @@ std::string plot (ImageType& image, int plot_axis)
ssize_t current_index = image.index (plot_axis);
image.index(plot_axis) = 0;

std::vector<value_type> plotslice (image.size(plot_axis));
std::vector<value_type> plotslice_finite (image.size(plot_axis));
std::vector<default_type> plotslice (image.size(plot_axis));
std::vector<default_type> plotslice_finite (image.size(plot_axis));
size_t k = 0;
for (auto l = Loop (plot_axis)(image); l; ++l) {
plotslice[k] = image.value();
plotslice[k] = get_scalar (image.value());
plotslice_finite[k] = plotslice[k];
++k;
}
value_type vmin = percentile(plotslice_finite, 0); // non-finite values removed
value_type vmax = percentile(plotslice_finite, 100);
default_type vmin = percentile(plotslice_finite, 0); // non-finite values removed
default_type vmax = percentile(plotslice_finite, 100);
if (vmax == vmin) {
vmin -= 1e-3;
vmax += 1e-3;
Expand Down Expand Up @@ -654,6 +690,7 @@ void show_help ()
+ key ("f", "show / hide crosshairs")
+ key ("r", "reset focus")
+ key ("i", "toggle between nearest (default) and linear interpolation")
+ key ("z", "toggle between real, imag, abs & arg display for complex numbers")
+ key ("left mouse & drag", "move focus")
+ key ("right mouse & drag", "adjust brightness / contrast")
+ key ("Esc", "reset brightness / contrast")
Expand Down Expand Up @@ -834,6 +871,7 @@ class CallBack : public EventLoop::CallBack
case 'r': focus[x_axis] = std::round (image.size(x_axis)/2); focus[x_axis] = std::round (image.size(x_axis)/2);
focus[slice_axis] = std::round (image.size(slice_axis)/2); break;
case 'i': interpolate = !interpolate; break;
case 'z': display_type++; if (display_type>3) display_type = 0; break;
case '+': zoom *= 1.1; std::cout << ClearScreen; break;
case '-': zoom /= 1.1; std::cout << ClearScreen; break;
case ' ':
Expand Down Expand Up @@ -888,7 +926,7 @@ class CallBack : public EventLoop::CallBack

void run ()
{
auto image = Image<value_type>::open (argument[0]);
auto image = ImageType::open (argument[0]);

size_t projection_axes[3] = {get_options("sagittal").size(), get_options("coronal").size(), get_options("axial").size()};
size_t psum = 0;
Expand All @@ -904,6 +942,8 @@ void run ()

int colourmap_ID = get_option_value ("colourmap", 0);

display_type = get_option_value ("complex", 0);

do_plot = get_options ("plot").size();
plot_axis = get_option_value ("plot", plot_axis);
if (plot_axis >= int(image.ndim()))
Expand Down