- args_info_clitkGammaIndex args_info;
-
- if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
- exit(1);
-
- if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
- std::cerr << "Specify either relative or absolute dose margin" << endl;
- exit(1);
- }
-
- if (args_info.isodose_number_arg <= 0) {
- std::cerr << "Specify a valid isodose number (>0)" << endl;
- exit(1);
- }
-
- bool verbose = args_info.verbose_flag;
-
- std::string reference_filename(args_info.reference_arg);
- std::string target_filename(args_info.target_arg);
- std::string gamma_filename(args_info.output_arg);
- float space_margin = args_info.spatial_margin_arg;
- float dose_rel_margin = args_info.relative_dose_margin_arg;
- float dose_margin = args_info.absolute_dose_margin_arg;
- bool use_dose_margin = args_info.absolute_dose_margin_given;
- unsigned int dose_size = args_info.isodose_number_arg;
-
- if (verbose) {
- cout << "reference_filename=" << reference_filename << endl;
- cout << "target_filename=" << target_filename << endl;
- cout << "gamma_filename=" << gamma_filename << endl;
- cout << "dose_size=" << dose_size << endl;
- cout << "space_margin=" << space_margin << endl;
- if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
- else cout << "dose_rel_margin=" << dose_rel_margin << endl;
- }
-
- // load images
- Reader::Pointer reference_reader = Reader::New();
- Reader::Pointer target_reader = Reader::New();
- {
- reference_reader->SetFileName(reference_filename);
- target_reader->SetFileName(target_filename);
- reference_reader->Update();
- target_reader->Update();
- }
-
- // intensity normalisation
- if (!use_dose_margin) {
- MinMaxer::PixelType reference_max = GetImageMaximum(reference_reader->GetOutput());
- //MinMaxer::PixelType target_max = GetImageMaximum(target_reader->GetOutput());
-
- dose_margin = reference_max*dose_rel_margin;
-
- if (verbose) cout << "dose_margin=" << dose_margin << endl;
- }
-
- // scale intensity
- Normalizer::Pointer reference_normalizer = Normalizer::New();
- Normalizer::Pointer target_normalizer = Normalizer::New();
- {
- reference_normalizer->SetShift(0);
- reference_normalizer->SetScale(1/dose_margin);
- reference_normalizer->SetInput(reference_reader->GetOutput());
- reference_normalizer->Update();
-
- target_normalizer->SetShift(0);
- target_normalizer->SetScale(1/dose_margin);
- target_normalizer->SetInput(target_reader->GetOutput());
- target_normalizer->Update();
-
- //cout << "scale=" << reference_normalizer->GetScale() << "/" << target_normalizer->GetScale() << endl;
- //cout << "shift=" << reference_normalizer->GetShift() << "/" << target_normalizer->GetShift() << endl;
- }
-
- // normalize space coordinates
- Scaler::Pointer reference_scaler = Scaler::New();
- Scaler::Pointer target_scaler = Scaler::New();
- {
- reference_scaler->SetInput(reference_normalizer->GetOutput());
- TuneScaler(reference_scaler,space_margin);
- reference_scaler->Update();
-
- target_scaler->SetInput(target_normalizer->GetOutput());
- TuneScaler(target_scaler,space_margin);
- target_scaler->Update();
-
- //SaveImage(reference_scaler->GetOutput(),"norm_reference.mhd");
- //SaveImage(reference_scaler->GetOutput(),"norm_target.mhd");
- }
-
- // compute hyper surface plane
- float reference_dose_scaled_max = GetImageMaximum(reference_scaler->GetOutput());
- float target_dose_scaled_max = GetImageMaximum(target_scaler->GetOutput());
- float dose_scaled_max = reference_dose_scaled_max > target_dose_scaled_max ? reference_dose_scaled_max : target_dose_scaled_max;
- ImageBin::Pointer image_bin = AllocateImageBin(reference_scaler->GetOutput(),target_scaler->GetOutput(),dose_size,dose_scaled_max);
- FillImageBin(image_bin,reference_scaler->GetOutput());
- //SaveImage(image_bin.GetPointer(),"surface.mhd");
-
- // compute distance map
- Mapper::Pointer mapper = Mapper::New();
- mapper->InsideIsPositiveOn();
- mapper->SquaredDistanceOff();
- mapper->UseImageSpacingOn();
- mapper->SetInput(image_bin);
- mapper->Update();
- //SaveImage(mapper->GetOutput(),"distance.mhd");
-
- // extract gamma index from distance map
- Image::Pointer image_gamma = AllocateImageGamma(target_normalizer->GetOutput());
- FillImageGamma(image_gamma,target_scaler->GetOutput(),mapper->GetOutput());
- SaveImage(image_gamma.GetPointer(),gamma_filename);
-
- if (verbose) ComputeGammaRatio(image_gamma);
-
- return 0;
-}
+ clitk::RegisterClitkFactories();
+
+ args_info_clitkGammaIndex args_info;
+
+ if (cmdline_parser_clitkGammaIndex(argc, argv, &args_info) != 0)
+ exit(1);
+
+ if (!args_info.absolute_dose_margin_given && !args_info.relative_dose_margin_given) {
+ std::cerr << "Specify either relative or absolute dose margin" << endl;
+ exit(1);
+ }
+
+ bool verbose = args_info.verbose_flag;
+
+ std::string reference_filename(args_info.reference_arg);
+ std::string target_filename(args_info.target_arg);
+ std::string gamma_filename(args_info.output_arg);
+ double space_margin = args_info.spatial_margin_arg;
+ double dose_rel_margin = args_info.relative_dose_margin_arg;
+ double dose_margin = args_info.absolute_dose_margin_arg;
+ bool use_dose_margin = args_info.absolute_dose_margin_given;
+
+ if (verbose) {
+ cout << "reference_filename=" << reference_filename << endl;
+ cout << "target_filename=" << target_filename << endl;
+ cout << "gamma_filename=" << gamma_filename << endl;
+ cout << "space_margin=" << space_margin << endl;
+ if (use_dose_margin) cout << "dose_margin=" << dose_margin << endl;
+ else cout << "dose_rel_margin=" << dose_rel_margin << endl;
+ }