1 /*=========================================================================
6 Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
7 l'Image). All rights reserved. See Doc/License.txt or
8 http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notices for more information.
14 =========================================================================*/
16 #include "clitkSignalMeanPositionFilter.h"
18 //---------------------------------------------------------------------
19 void clitk::SignalMeanPositionFilter::SetParameters(args_info_clitkSignalMeanPositionTracking & args) {
22 mOutputFilenameIsSet = false;
23 mOutputResidualFilenameIsSet = false;
24 mInput.Read(args_info.input_arg);
25 mAugmentationDelay = args_info.delay_arg;
26 mIsAdaptiveMethod = false;
28 if (args_info.eta_given) {
29 mEta = args_info.eta_arg;
32 mMaxIteration = args_info.iter_arg;
33 if (args_info.output_given) {
34 mOutputFilenameIsSet = true;
35 mOutputFilename = args_info.output_arg;
37 if (args_info.residual_given) {
38 mOutputResidualFilenameIsSet = true;
39 mOutputResidualFilename = args_info.residual_arg;
41 if (args_info.augmented_given) {
42 mOutputAugmentedFilenameIsSet = true;
43 mOutputAugmentedFilename = args_info.augmented_arg;
45 mVerbose = args_info.verbose_flag;
46 mVerboseIteration = args_info.verbose_iteration_flag;
47 if (args_info.L_given) {
48 mIsAdaptiveMethod = true;
49 mWindowLength = args_info.L_arg;
52 //---------------------------------------------------------------------
55 //---------------------------------------------------------------------
56 void clitk::SignalMeanPositionFilter::Update() {
60 // clitk::Signal temp;
61 // temp.resize(mInput.size()*e);
62 // for(unsigned int i=0; i<temp.size(); i++) {
63 // int index = lrint(floor((double)i/e));
64 // double res = (double)i/e - index;
67 // temp[i] = mInput[index] + res*(mInput[index+1]-mInput[index])/e;
69 // mInput.resize(temp.size());
70 for (unsigned int i=0; i<mInput.size(); i++) mInput[i] = mInput[i]+(double)i/(double)mInput.size();
72 for (unsigned int i=0; i<mInput.size(); i++) {
73 mInput[i] = (0.8+((double)rand()/ (double)RAND_MAX)*0.4) * mInput[i];
78 openFileForWriting(os, "input.dat");
79 for(unsigned int i=0; i<mInput.size(); i++) os << mInput[i] << std::endl;
84 // Compute augmented space
86 std::cout << "Computing augmented space with delay = " << mAugmentationDelay
87 << ", number of samples is " << mInput.size()-mAugmentationDelay
90 ComputeAugmentedSpace(mInput, mAugmentedInputX, mAugmentedInputY, mAugmentationDelay);
92 // Output augmented state
93 if (mOutputAugmentedFilenameIsSet) {
95 openFileForWriting(os, mOutputAugmentedFilename);
96 for(unsigned int i=0; i<(unsigned int)mWindowLength; i++) {
97 os << mAugmentedInputX[i] << " " << mAugmentedInputY[i] << std::endl;
102 // Initial starting ellipse
103 Ellipse An; // starting point
104 if (!mEtaIsSet) mEta = -1;
105 An.InitialiseEllipseFitting(mEta, mWindowLength, mAugmentedInputX, mAugmentedInputY);
108 std::cout << "Eta is " << mEta << std::endl;
110 Ellipse AnInitial(An);
112 std::cout << "Set initial ellipse to c=" << An.ComputeCenter()
113 << " and axes=" << An.ComputeSemiAxeLengths() << std::endl;
117 if (mIsAdaptiveMethod) {
118 AdaptiveFitEllipse(An);
125 std::cout << "Result is " << An << std::endl;
126 std::cout << "Center is " << An.ComputeCenter() << std::endl;
127 std::cout << "SemiAxis are " << An.ComputeSemiAxeLengths() << std::endl;
128 std::cout << "Angle is " << rad2deg(An.ComputeAngleInRad()) << " deg" << std::endl;
130 if (mOutputFilenameIsSet) {
132 openFileForWriting(os, mOutputFilename);
133 os << AnInitial.ComputeCenter()[0] << " " << AnInitial.ComputeCenter()[1] << " "
134 << AnInitial.ComputeSemiAxeLengths()[0] << " " << AnInitial.ComputeSemiAxeLengths()[1] << " "
135 << AnInitial.ComputeAngleInRad();
136 for(int i=0; i<6; i++) os << AnInitial[i] << " " ;
139 os << An.ComputeCenter()[0] << " " << An.ComputeCenter()[1] << " "
140 << An.ComputeSemiAxeLengths()[0] << " " << An.ComputeSemiAxeLengths()[1] << " "
141 << An.ComputeAngleInRad();
142 for(int i=0; i<6; i++) os << An[i] << " ";
146 openFileForWriting(os, "centers.dat");
147 for(unsigned int i=0; i<mCenters.size(); i++) {
148 os << mCenters[i][0] << " " << mCenters[i][1] << std::endl;
154 if (mOutputResidualFilenameIsSet) {
156 openFileForWriting(os, mOutputResidualFilename);
157 for(unsigned int i=0; i<mCurrentResidual.size(); i++) {
159 os << mCurrentResidual[i] << std::endl;
164 //---------------------------------------------------------------------
167 //---------------------------------------------------------------------
168 void clitk::SignalMeanPositionFilter::ComputeAugmentedSpace(const clitk::Signal & input,
169 clitk::Signal & outputX,
170 clitk::Signal & outputY,
171 unsigned int delay) {
172 if (input.size() <= delay) {
173 std::cerr << "Error in signal length is " << input.size()
174 << " while delay is " << delay << " : too short. Abort." << std::endl;
177 outputX.resize(input.size()-delay);
178 outputY.resize(input.size()-delay);
179 for(unsigned int i=0; i<outputX.size(); i++) {
180 outputX[i] = input[i+delay];
181 outputY[i] = input[i];
185 //---------------------------------------------------------------------
188 //---------------------------------------------------------------------
189 void clitk::SignalMeanPositionFilter::FitEllipse(clitk::Ellipse & An) {
191 mCurrentResidual.clear();
192 while (n<mMaxIteration) {
193 double r = An.EllipseFittingNextIteration();
194 mCurrentResidual.push_back(r);
196 if (mVerboseIteration) {
197 std::cout.precision(3);
198 std::cout << "A(" << n << ")=" << An
199 << " c=" << An.ComputeCenter()
200 << " ax=" << An.ComputeSemiAxeLengths()
201 << " theta=" << rad2deg(An.ComputeAngleInRad());
202 std::cout.precision(10);
203 std::cout << " R=" << r << std::endl;
207 //---------------------------------------------------------------------
210 //---------------------------------------------------------------------
211 void clitk::SignalMeanPositionFilter::AdaptiveFitEllipse(clitk::Ellipse & An) {
212 for(unsigned int t=0; t<(unsigned int)args_info.t_arg; t++) {
214 An.UpdateSMatrix(t, mWindowLength, mAugmentedInputX, mAugmentedInputY);
215 mCenters.push_back(An.ComputeCenter());
219 //---------------------------------------------------------------------