2 ***********************************************************
3 $SNARK_Header: S N A R K 1 4 - A PICTURE RECONSTRUCTION PROGRAM $
4 $HeadURL: svn://dig.cs.gc.cuny.edu/snark/trunk/src/snark/mart.cpp $
5 $LastChangedRevision: 80 $
6 $Date: 2014-07-01 21:01:54 -0400 (Tue, 01 Jul 2014) $
8 ***********************************************************
10 IMPLEMENTATIONS OF VARIOUS VERSIONS OF MULTIPLICATIVE ART ALGORITHM
11 PROGRAMMED BY A. V. LAKSHMINARAYANAN AND A. LENT
12 SUNY AT BUFFALO, MAY, 1977; REVISED JAN 78
13 IT IS ASSUMED THAT PROJECTION DATA ARE INTRINSICALLY NONNEGATIVE
14 NONNEGATIVE RECONSTRUCTIONS ARE SOUGHT
16 DESCRIPTION OF PARAMETERS
18 INPUT CONSISTS OF THE KEYWORD 'METHOD' FOLLOWED BY ONE OF FOUR
19 OPTIONS 'GBH', 'LENT', 'LAK1', 'LAK2', FOLLOWED BY AN INTEGER
20 MODIFIER KOUNT, TWO FLOATING POINT MODIFIERS RELAX AND DATZRO,
21 OPTIONALLY FOLLOWED BY WORD MODIFIERS 'NORMAL', 'ENTROPY'
22 THE KEYWORD 'METHOD' MUST BE PRESENT; ITS MODIFIERS HAVE THE
25 'GBH', USE MULTIPLICATIVE ART ALGORITHM OF GORDON AND HERMAN
26 INT. REV. CYTOL. 38, 1974
27 START WITH PICTURE UNIFORMLY AVEDEN
28 RELAXATION PARAMETER ADDED TO LITERATURE VERSION OF
30 'LENT', USE MULTIPLICATIVE ART ALGORITHM OF LENT
31 PROC. CONF ON IMAGE ANALYSIS AND EVALUATION, JULY,1976, TORON
32 SOCIETY OF PHOTOGRAPHIC SCIENTISTS AND ENGINEERS
33 START WITH PICTURE AVEDEN EXCEPT FOR PIXELS IN RAYS WITH ZERO
35 THIS STARTING POINT IS IMMEDIATELY CHANGED BY APPLICATION
36 OF NORMALIZATION CONSTRAINT
37 NORMALIZE COEFFICIENT MATRIX BY ROWS
39 'LAK1', USE (APPROX) MULTIPLICATIVE ART ALGORITHM OF LAKSHMINARAYANAN
40 LINEAR APPROXIMATION TO THE NONINTEGRAL POWER OF MULT ART.
41 FIRST TWO TERMS OF THE BINOMIAL EXPANSION USED.
42 USE SAME STARTING POINT AS UNDER 'LENT'
43 NORMALIZE COEFFICIENTS BY MAX ELEMENT OF PROJECTION MATRIX
45 'LAK2', USE (APPROX) MULT ART OF LAKSHMINARAYANAN
46 QUADRATIC APPROXIMATION TO THE NONINTEGRAL POWER OF MULT ART.
47 FIRST THREE TERMS OF THE BINOMIAL EXPANSION USED.
48 SAME STARTING POINT AND COEFFICIENT NORMALIZATION AS 'LAK1'
50 IF KOUNT .NE. 0, USE ONLY THE FIRST KOUNT RAYS OF THE PICK SEQUENCE
51 THE EFFECT OF KOUNT IS DETERMINED BY THE SNARK SELECT CARD
52 DEFAULT...IF KOUNT .LE. 0, KOUNT RESET TO COMPLETE SET OF PROJ DATA
54 RELAX IS THE UNDERRELXATION PARAMETER.
55 VALUES LESS THAN 1.0 ARE SUGGESTED FOR NOISY PROJECTION DATA.
56 VALUES BIGGER THAN 1.0 ARE NOT PROHIBITED BY THEORY AND MAY SPEED
58 THE PRECEDING REMARKS ARE NOT MEANINGFUL FOR 'GBH'
59 FOR 'GBH', RELAX MUST BE CONSIDERED TO HAVE DIMENSIONS
61 DATZRO..ANY PROJECTION DATA WITH VALUE .LE. THAN DATZRO WILL BE
62 CONSIDERED ZERO BY THE PROGRAM.
64 'NORMAL' , IF PRESENT, NORMALIZE RECON TO DENSITY AVEDEN
65 AT END OF EACH ITERATION
67 'ENTROPY' ,IF PRESENT, CALCULATE AND PRINT ENTROPY FUNCTIONAL
69 LET T = SIGMA RECON(I) OVER ALL PIXELS
71 ENTROPY = - SIGMA Y(I)*ALOG(Y(I))
73 EXAMPLE INPUT SEQUENCES FOLLOW
74 METHOD GBH 0 RELAX = 1.0 IGNORE RAYSUMS LT 0.0
75 METHOD LAK1 KOUNT = 0 0.8 0.0 ENTROPY
95 #include "projfile.h" //for debugg purpose
99 BOOLEAN mart_class::Run(REAL* recon, INTEGER* list, REAL* weight, INTEGER iter)
101 static const INTEGER hmeth = CHAR2INT('m', 'e', 't', 'h');
102 static const INTEGER hnorm = CHAR2INT('n', 'o', 'r', 'm');
103 static const INTEGER hentr = CHAR2INT('e', 'n', 't', 'r');
105 static const INTEGER gbh = CHAR2INT('g', 'b', 'h', ' ');
106 static const INTEGER lent = CHAR2INT('l', 'e', 'n', 't');
107 static const INTEGER lak1 = CHAR2INT('l', 'a', 'k', '1');
108 static const INTEGER lak2 = CHAR2INT('l', 'a', 'k', '2');
114 static REAL datzro; //set static
115 static INTEGER lnorm;
116 static INTEGER lfnal;
127 static REAL aijmax; // set static
144 static INTEGER kount;
145 static REAL relax; //set static
152 // num_rays = GeoPar.nrays;
158 // num_rays = GeoPar.usrays;
170 word = InFile.getwrd(TRUE, &eol, exp0, 1);
175 "\n **** error - keyword method is missing***");
180 { gbh, lent, lak1, lak2 };
182 word = InFile.getwrd(FALSE, &eol, exp1, 4);
186 fprintf(output, "\n **** error - one and only one of the following methods must be specified : gbh, lent, lak1, or lak2***");
193 fprintf(output, "\n multiplicative art version %s",
196 else if (word == lent)
199 fprintf(output, "\n multiplicative art version %s",
202 else if (word == lak1)
205 fprintf(output, "\n multiplicative art version %s",
211 fprintf(output, "\n multiplicative art version %s",
215 // GET KOUNT, RELAXATION PARAMETER, RAYSUM LIMIT, NORMAL, ENTROPY
216 kount = InFile.getint(FALSE, &eol);
218 relax = InFile.getnum(FALSE, &eol);
220 datzro = InFile.getnum(FALSE, &eol);
225 "\n **** some data is missing: either kount, relax, or data_zero***");
233 word = InFile.getwrd(FALSE, &eol, exp3, 2);
237 // finished reading input
247 word = InFile.getwrd(FALSE, &eol, exp3a, 1);
259 //*** CHECKING INPUTS AND SETTING DEFAULTS
262 if (RaySel.useray) // select user
263 kount = GeoPar.prjnum * GeoPar.usrays; // change nrays to usrays. hstau 7/2003
266 kount = GeoPar.prjnum * GeoPar.snrays;
267 } // set value to kount for select user or select snark. Wei 11/2004
271 if ((nalg == 3) && (relax > 1.0))
273 if (datzro < Consts.zero)
274 datzro = Consts.zero;
276 fprintf(output, "\n no of rays used in each iteration %10i",
278 fprintf(output, "\n underrelaxation factor %10.5f", relax);
280 "\n all projection data with values less than %20.10e are ignored",
285 fprintf(output, "\n reconstruction is normalized");
289 fprintf(output, "\n reconstruction not normalized");
294 fprintf(output, "\n entropy functional is calculated");
298 fprintf(output, "\n entropy functional not calculated");
301 if (GeoPar.aveden <= Consts.zero)
304 "\n ***aveden too small. probably incorrect use of mart, no execution***");
308 for (i = 0; i < area; i++)
310 recon[i] = GeoPar.aveden;
313 //*** NO CHOICE IN INITIAL GUESS; EXALG OPTIONS DONT OPERATE
317 // WILL ZERO OUT ANY PIXEL IN A RAY WITH ZERO PROJ DATA
318 for (np = 0; np < GeoPar.prjnum; np++)
320 for (nr = 0; nr < GeoPar.nrays; nr++)
322 raysum = Anglst.prdta(np, nr);
323 if (raysum <= datzro)
329 ray(np, nr, list, weight, &numb, &snorm);
333 Blob.bray(np, nr, list, weight, &numb, &snorm);
340 wray(np, nr, list, weight, &numb, &snorm);
344 Blob.bwray(np, nr, list, weight, &numb, &snorm);
349 for (j = 0; j < numb; j++)
362 for (i = 0; i < area; i++)
369 if (i % 2 == Blob.pr)
385 "\n ***warning- all projection data near zero. zero reconstruction.***");
391 fact = area * GeoPar.aveden / isum;
396 for (i = 0; i < area; i++)
402 // LOCATE BIGGEST ELEMENT IN PROJECTION MATRIX FOR LAK1 AND LAK2 ONLY
403 // THIS STRATEGY IS CONSERVATIVE IN THAT COEFFICIENTS OF PIXELS
404 // WHOSE VALUES ARE KNOWN TO BE ZERO ARE INCLUDED IN THE SEARCH
408 for (np = 0; np < GeoPar.prjnum; np++)
410 for (nr = 0; nr < GeoPar.nrays; nr++)
416 ray(np, nr, list, weight, &numb, &snorm);
420 Blob.bray(np, nr, list, weight, &numb, &snorm);
427 wray(np, nr, list, weight, &numb, &snorm);
431 Blob.bwray(np, nr, list, weight, &numb, &snorm);
436 for (j = 0; j < numb; j++)
447 // BASIC UPDATE PART OF PROGRAM
455 // GORDON,BENDER,HERMAN MULT ART
457 for (k = 0; k < kount; k++)
461 raysum = Anglst.prdta(np, nr);
464 trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
469 trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
470 &snorm, GeoPar.line, FALSE);
474 if (raysum <= datzro)
476 for (j = 0; j < numb; j++)
485 if ((trysum <= Consts.zero) && (trace >= 1))
488 "\n zero pseudo proj dat %5i %5i %10.5f %10.5f",
489 np, nr, raysum, trysum);
492 if (trysum > Consts.zero)
494 corfac = relax * (REAL) log(raysum / trysum);
495 for (j = 0; j < numb; j++)
498 recon[l] *= (REAL) exp(corfac * weight[j]);
512 for (k = 0; k < kount; k++)
515 raysum = Anglst.prdta(np, nr);
520 trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
525 trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
526 &snorm, GeoPar.line, FALSE);
532 delu = (REAL) log(raysum / trysum);
534 // FIND ROW SCALE FACTOR
537 for (j = 0; j < numb; j++)
539 if (weight[j] > ajmax)
543 corr = relax * delu / ajmax;
545 for (j = 0; j < numb; j++)
548 recon[l] *= (REAL) exp(corr * weight[j]);
559 // LAK APPROX MULT ART - LINEAR APPROXIMATION
562 roaij = relax / aijmax;
564 for (k = 0; k < kount; k++)
567 raysum = Anglst.prdta(np, nr);
572 trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
577 trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
578 &snorm, GeoPar.line, FALSE);
580 if (!((numb == 0) || (trysum <= Consts.zero)))
582 reldif = (raysum - trysum) / MAX0(raysum, trysum);
585 for (j = 0; j < numb; j++)
588 aij = weight[j] * roaij;
589 recon[l] *= ((REAL) 1.0 + reldif * aij);
597 // LAK APPROX MULT ART - QUADRATIC APPROXIMATION
601 roaij = relax / aijmax;
603 for (k = 0; k < kount; k++)
606 raysum = Anglst.prdta(np, nr);
611 trysum = pseudo(recon, np, nr, list, weight, &numb, &snorm,
616 trysum = Blob.bpseudo(recon, np, nr, list, weight, &numb,
617 &snorm, GeoPar.line, FALSE);
619 if (!((numb == 0) || (trysum <= Consts.zero)))
621 reldif = (raysum - trysum) / MAX0(raysum, trysum);
624 hafrel = (REAL) 0.5 * reldif;
625 hafabs = (REAL) fabs(hafrel);
626 for (j = 0; j < numb; j++)
629 aij = weight[j] * roaij;
632 * ((REAL) 1.0 + hafrel * aij + hafabs);
645 if (!(lnorm || lfnal))
650 for (i = 0; i < area; i++)
654 if (i % 2 == Blob.pr)
665 if (tot > Consts.zero)
669 fact = GeoPar.aveden * area / tot; // aveden in both basis concide
674 for (i = 0; i < area; i++)
679 tot = GeoPar.aveden * area;
690 for (i = 0; i < area; i++)
694 rot = recon[i] / tot;
695 entrop -= rot * (REAL) log(rot); //remove the minus sign here, wei
700 "\n iter %5i entropy %20.10e",
706 "\n ***total density %20.10e is too small, no normalization or entropy calculation***",