2 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
3 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/superior.cpp $
4 $LastChangedRevision: 158 $
5 $Date: 2014-07-30 16:59:01 -0400 (Wed, 30 Jul 2014) $
43 sec_cri_class* scClasses[] =
58 nav_class* navClasses[] =
67 class superior_class SuperSet =
79 void initSuperiorization()
85 // enable superiorization
86 SuperSet.superiorizationEnabled = TRUE;
87 cout << "\n Superiorization is enabled";
90 word = InFile.getint(FALSE, &eol);
94 cout << "\n **** N must be greater than zero";
95 cout << "\n **** program aborted\n";
99 cout << "\n N = " << SuperSet.N;
102 number = InFile.getnum(FALSE, &eol);
106 cout << "\n **** a must be greater than zero";
107 cout << "\n **** program aborted\n";
112 cout << "\n **** a must be less than one";
113 cout << "\n **** program aborted\n";
117 cout << "\n a = " << SuperSet.a;
120 number = InFile.getnum(FALSE, &eol);
124 cout << "\n **** b must be greater than zero";
125 cout << "\n **** program aborted\n";
129 cout << "\n b = " << SuperSet.b;
131 // read secondary criterion
132 static const INTEGER criterion[] =
134 CHAR2INT('t', 'v', 'a', 'r'),
135 CHAR2INT('s', 'm', 'o', 'o'),
136 CHAR2INT('s', 'c', 'r', '3'),
137 CHAR2INT('s', 'c', 'r', '4'),
138 CHAR2INT('s', 'c', 'r', '5')
141 word = InFile.getwrd(FALSE, &eol, criterion, 5);
143 for (INTEGER i = 0; i < 5; i++)
145 if (word == criterion[i])
147 cout << "\n secondary criterion: " << int2str(criterion[i]);
149 SuperSet.secondaryCriterion = scClasses[i];
150 SuperSet.secondaryCriterion->Init();
152 // set the default non-ascending vector for the selected secondary criterion
153 SuperSet.nonascendingVector = navClasses[i];
154 SuperSet.nonascendingVector->Init();
160 if (SuperSet.secondaryCriterion == NULL)
162 cout << "\n **** unknown secondary criterion ", int2str(word);
163 cout << "\n **** program aborted\n";
167 // set default values of the optional commands
168 SuperSet.report = false;
169 SuperSet.alternativel=0;
170 SuperSet.positivity = false;
172 // read remaining optional commands: posi, altl (alternative l) or rprt
173 static const INTEGER optionsf[] =
175 CHAR2INT('a', 't', 'l', '1'),
176 CHAR2INT('a', 't', 'l', '2'),
177 CHAR2INT('r', 'p', 'r', 't'),
178 CHAR2INT('p', 'o', 's', 'i')
181 word = InFile.getwrd(FALSE, &eol, optionsf, 4);
183 for (INTEGER i = 0; i < 4; i++)
185 if (word == optionsf[i])
189 SuperSet.alternativel=1;
190 cout << "\n alternative l handling 1 is enabled";
194 SuperSet.alternativel=2;
195 cout << "\n alternative l handling 2 is enabled";
199 SuperSet.report = true;
200 cout << "\n reporting is enabled";
201 cout << "\n reporting file: RPRTsuperiorization";
203 SuperSet.skips = InFile.getnum(FALSE, &eol);
204 if (SuperSet.skips>0)
206 cout << "\n report skip factor " << SuperSet.skips;
210 cout << "\n reporting on every iteration";
215 SuperSet.positivity=true;
216 cout << "\n positivity constraint is enabled";
221 // read remaining optional commands: altl (alternative l) or rprt
222 static const INTEGER options[] =
224 CHAR2INT('a', 't', 'l', '1'),
225 CHAR2INT('a', 't', 'l', '2'),
226 CHAR2INT('r', 'p', 'r', 't')
229 word = InFile.getwrd(FALSE, &eol, options, 3);
231 for (INTEGER i = 0; i < 3; i++)
233 if (word == options[i])
237 SuperSet.alternativel=1;
238 cout << "\n alternative l handling 1 is enabled";
242 SuperSet.alternativel=2;
243 cout << "\n alternative l handling 2 is enabled";
247 SuperSet.report = true;
248 cout << "\n reporting is enabled";
249 cout << "\n reporting file: RPRTsuperiorization";
251 SuperSet.skips = InFile.getnum(FALSE, &eol);
252 if (SuperSet.skips>0)
254 cout << "\n report skip factor " << SuperSet.skips;
258 cout << "\n reporting on every iteration";
264 // read last possible command (rprt)
265 INTEGER rprt = CHAR2INT('r','p','r','t');
266 if (InFile.getwrd(FALSE, &eol, &rprt, 1) == rprt)
268 SuperSet.report = true;
269 cout << "\n reporting is enabled";
270 cout << "\n reporting file: RPRTsuperiorization";
272 SuperSet.skips = InFile.getnum(FALSE, &eol);
273 if (SuperSet.skips>0)
275 cout << "\n report skip factor " << SuperSet.skips;
279 cout << "\n reporting on every iteration";
287 BOOLEAN executeSuperiorization(REAL* recon, INTEGER* list, REAL* weight, INTEGER count, alg_class* algorithm)
289 double min = pow(10,-15);
294 pixelcount = GeoPar.area;
298 cout << "\nError: Superiorization does currently not support blobs, skipping Superiorization.";
300 // execute one iteration of the selected algorithm
301 return algorithm->Run(recon, list, weight, count);
304 // user defined superiorization parameters
305 INTEGER N = SuperSet.N;
308 sec_cri_class* secCrit = SuperSet.secondaryCriterion;
309 nav_class* naV = SuperSet.nonascendingVector;
311 // xkn, vkn and z are declared here so that memory is only allocated once
312 REAL* xkn = new REAL[pixelcount];
313 REAL* vkn = new REAL[pixelcount];
314 REAL* z = new REAL[pixelcount];
316 // things that should only be done on the first iteration
318 // taken care of by loop in exalg
320 // not necessary (accessing recon directly)
323 if (SuperSet.alternativel==1)
327 else if (SuperSet.alternativel==2 && count>1)
329 l = (int)round(Rand()*(SuperSet.l-count)) + count - 1;
341 // taken care of by loop in exalg
345 memcpy(xkn, recon, pixelcount * sizeof(REAL));
347 REAL phixk=secCrit->Run(recon); // only needs to be processed once
350 // "set vkn to be a nonascending vector for phi at xkn"
360 // "set betakn = mul"
361 REAL betakn = b*pow(a, l);
362 if (betakn<min) betakn=0;
363 // "set z = xkn + betakn*vkn"
364 for (INTEGER i = 0; i < pixelcount; i++) z[i] = xkn[i] + (betakn * vkn[i]);
365 // "if z e delta and phi(z) <= phi(xk) then"
366 BOOLEAN inDelta = true;
367 if (SuperSet.positivity)
369 for (INTEGER i = 0; i < pixelcount; i++) if (z[i] < 0) inDelta=false;
371 if (inDelta && secCrit->Run(z) <= phixk)
376 memcpy(xkn, z, pixelcount * sizeof(REAL));
377 // "set loop = false"
382 // "set xk+1 = B*xkN"
383 memcpy(recon, xkn, pixelcount * sizeof(REAL));
389 // save value of l for next iteration (checking if this is necessary will take as long as saving it)
393 if (SuperSet.report && (count==1 || SuperSet.skips==0 || (count%(INTEGER)SuperSet.skips)==0))
395 REAL preAlgo = secCrit->Run(recon);
396 BOOLEAN error = algorithm->Run(recon, list, weight, count);
397 cout << "\n value of l: " << l;
398 cout << "\n value of phi before algorithm operator: " << preAlgo;
399 REAL postAlgo = secCrit->Run(recon);
400 cout << "\n value of phi after algorithm operator: " << postAlgo;
405 RPRTfile = fopen("RPRTsuperiorization", "w");
406 fprintf(RPRTfile,"superiorization reporting output");
407 fprintf(RPRTfile,"\n iter l phi pre-algorithm phi post-algorithm");
408 fprintf(RPRTfile,"\n %5i %10i %20.9f %20.9f", count, l, preAlgo, postAlgo);
412 RPRTfile = fopen("RPRTsuperiorization", "a");
413 fprintf(RPRTfile,"\n %5i %10i %20.9f %20.9f", count, l, preAlgo, postAlgo);
420 // execute one iteration of the selected algorithm
421 return algorithm->Run(recon, list, weight, count);
424 // taken care of by loop in exalg