gpfa3f.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

/*<       if (ja.lt.istart) ja = ja + ninc >*/
                        if (ja < istart) {
                            ja += ninc;
                        }
/*<   415 continue >*/
/* L415: */
                    }
/*<   420 continue >*/
/* L420: */
                }
/*<   430 continue >*/
/* L430: */
            }
/* -----( end of double loop for k=0 ) */

/*  finished if last pass */
/*  --------------------- */
/*<       if (ipass.eq.m) go to 490 >*/
            if (ipass == m) {
                goto L490;
            }

/*<       kk = 2*la >*/
            kk = la << 1;

/*     loop on nonzero k */
/*     ----------------- */
/*<       do 470 k = ink , jstep-ink , ink >*/
            i__4 = jstep - ink;
            i__3 = ink;
            for (k = ink; i__3 < 0 ? k >= i__4 : k <= i__4; k += i__3) {
/*<       co1 = trigs(kk+1) >*/
                co1 = trigs[kk + 1];
/*<       si1 = s*trigs(kk+2) >*/
                si1 = s * trigs[kk + 2];
/*<       co2 = trigs(2*kk+1) >*/
                co2 = trigs[(kk << 1) + 1];
/*<       si2 = s*trigs(2*kk+2) >*/
                si2 = s * trigs[(kk << 1) + 2];

/*  double loop along first transform in block */
/*  ------------------------------------------ */
/*<       do 460 ll = k , (la-1)*ink , 3*jstep >*/
                i__5 = (la - 1) * ink;
                i__6 = jstep * 3;
                for (ll = k; i__6 < 0 ? ll >= i__5 : ll <= i__5; ll += i__6) {

/*<       do 450 jjj = ll , (n-1)*inc , 3*la*ink >*/
                    i__7 = (*n - 1) * *inc;
                    i__8 = la * 3 * ink;
                    for (jjj = ll; i__8 < 0 ? jjj >= i__7 : jjj <= i__7; jjj 
                            += i__8) {
/*<       ja = istart + jjj >*/
                        ja = istart + jjj;

/*  "transverse" loop */
/*  ----------------- */
/*<       do 445 nu = 1 , inq >*/
                        i__9 = inq;
                        for (nu = 1; nu <= i__9; ++nu) {
/*<       jb = ja + jstepl >*/
                            jb = ja + jstepl;
/*<       if (jb.lt.istart) jb = jb + ninc >*/
                            if (jb < istart) {
                                jb += ninc;
                            }
/*<       jc = jb + jstepl >*/
                            jc = jb + jstepl;
/*<       if (jc.lt.istart) jc = jc + ninc >*/
                            if (jc < istart) {
                                jc += ninc;
                            }
/*<       jd = ja + laincl >*/
                            jd = ja + laincl;
/*<       if (jd.lt.istart) jd = jd + ninc >*/
                            if (jd < istart) {
                                jd += ninc;
                            }
/*<       je = jd + jstepl >*/
                            je = jd + jstepl;
/*<       if (je.lt.istart) je = je + ninc >*/
                            if (je < istart) {
                                je += ninc;
                            }
/*<       jf = je + jstepl >*/
                            jf = je + jstepl;
/*<       if (jf.lt.istart) jf = jf + ninc >*/
                            if (jf < istart) {
                                jf += ninc;
                            }
/*<       jg = jd + laincl >*/
                            jg = jd + laincl;
/*<       if (jg.lt.istart) jg = jg + ninc >*/
                            if (jg < istart) {
                                jg += ninc;
                            }
/*<       jh = jg + jstepl >*/
                            jh = jg + jstepl;
/*<       if (jh.lt.istart) jh = jh + ninc >*/
                            if (jh < istart) {
                                jh += ninc;
                            }
/*<       ji = jh + jstepl >*/
                            ji = jh + jstepl;
/*<       if (ji.lt.istart) ji = ji + ninc >*/
                            if (ji < istart) {
                                ji += ninc;
                            }
/*<       j = 0 >*/
                            j = 0;

/*  loop across transforms */
/*  ---------------------- */
/* dir$ ivdep, shortloop */
/*<       do 440 l = 1 , nvex >*/
                            i__10 = nvex;
                            for (l = 1; l <= i__10; ++l) {
/*<       ajb = a(jb+j) >*/
                                ajb = a[jb + j];
/*<       ajc = a(jc+j) >*/
                                ajc = a[jc + j];
/*<       t1 = ajb + ajc >*/
                                t1 = ajb + ajc;
/*<       aja = a(ja+j) >*/
                                aja = a[ja + j];
/*<       t2 = aja - 0.5 * t1 >*/
                                t2 = aja - t1 * (float).5;
/*<       t3 = c1 * ( ajb - ajc ) >*/
                                t3 = c1 * (ajb - ajc);
/*<       ajd = a(jd+j) >*/
                                ajd = a[jd + j];
/*<       ajb =  ajd >*/
                                ajb = ajd;
/*<       bjb = b(jb+j) >*/
                                bjb = b[jb + j];
/*<       bjc = b(jc+j) >*/
                                bjc = b[jc + j];
/*<       u1 = bjb + bjc >*/
                                u1 = bjb + bjc;
/*<       bja = b(ja+j) >*/
                                bja = b[ja + j];
/*<       u2 = bja - 0.5 * u1 >*/
                                u2 = bja - u1 * (float).5;
/*<       u3 = c1 * ( bjb - bjc ) >*/
                                u3 = c1 * (bjb - bjc);
/*<       bjd = b(jd+j) >*/
                                bjd = b[jd + j];
/*<       bjb =  bjd >*/
                                bjb = bjd;
/*<       a(ja+j) = aja + t1 >*/
                                a[ja + j] = aja + t1;
/*<       b(ja+j) = bja + u1 >*/
                                b[ja + j] = bja + u1;
/*<       a(jd+j) = co1*(t2-u3) - si1*(u2+t3) >*/
                                a[jd + j] = co1 * (t2 - u3) - si1 * (u2 + t3);
/*<       b(jd+j) = si1*(t2-u3) + co1*(u2+t3) >*/
                                b[jd + j] = si1 * (t2 - u3) + co1 * (u2 + t3);
/*<       ajc =  co2*(t2+u3) - si2*(u2-t3) >*/
                                ajc = co2 * (t2 + u3) - si2 * (u2 - t3);
/*<       bjc =  si2*(t2+u3) + co2*(u2-t3) >*/
                                bjc = si2 * (t2 + u3) + co2 * (u2 - t3);
/* ---------------------- */
/*<       aje = a(je+j) >*/
                                aje = a[je + j];
/*<       ajf = a(jf+j) >*/
                                ajf = a[jf + j];
/*<       t1 = aje + ajf >*/
                                t1 = aje + ajf;
/*<       t2 = ajb - 0.5 * t1 >*/
                                t2 = ajb - t1 * (float).5;
/*<       t3 = c1 * ( aje - ajf ) >*/
                                t3 = c1 * (aje - ajf);
/*<       ajh = a(jh+j) >*/
                                ajh = a[jh + j];
/*<       ajf =  ajh >*/
                                ajf = ajh;
/*<       bje = b(je+j) >*/
                                bje = b[je + j];
/*<       bjf = b(jf+j) >*/
                                bjf = b[jf + j];
/*<       u1 = bje + bjf >*/
                                u1 = bje + bjf;
/*<       u2 = bjb - 0.5 * u1 >*/
                                u2 = bjb - u1 * (float).5;
/*<       u3 = c1 * ( bje - bjf ) >*/
                                u3 = c1 * (bje - bjf);
/*<       bjh = b(jh+j) >*/
                                bjh = b[jh + j];
/*<       bjf =  bjh >*/
                                bjf = bjh;
/*<       a(jb+j) = ajb + t1 >*/
                                a[jb + j] = ajb + t1;
/*<       b(jb+j) = bjb + u1 >*/
                                b[jb + j] = bjb + u1;
/*<       a(je+j) = co1*(t2-u3) - si1*(u2+t3) >*/
                                a[je + j] = co1 * (t2 - u3) - si1 * (u2 + t3);
/*<       b(je+j) = si1*(t2-u3) + co1*(u2+t3) >*/
                                b[je + j] = si1 * (t2 - u3) + co1 * (u2 + t3);
/*<       a(jh+j) = co2*(t2+u3) - si2*(u2-t3) >*/
                                a[jh + j] = co2 * (t2 + u3) - si2 * (u2 - t3);
/*<       b(jh+j) = si2*(t2+u3) + co2*(u2-t3) >*/
                                b[jh + j] = si2 * (t2 + u3) + co2 * (u2 - t3);
/* ---------------------- */
/*<       aji = a(ji+j) >*/
                                aji = a[ji + j];
/*<       t1 = ajf + aji >*/
                                t1 = ajf + aji;
/*<       ajg = a(jg+j) >*/
                                ajg = a[jg + j];
/*<       t2 = ajg - 0.5 * t1 >*/
                                t2 = ajg - t1 * (float).5;
/*<       t3 = c1 * ( ajf - aji ) >*/
                                t3 = c1 * (ajf - aji);
/*<       t1 = ajg + t1 >*/
                                t1 = ajg + t1;
/*<       a(jg+j) = ajc >*/
                                a[jg + j] = ajc;
/*<       bji = b(ji+j) >*/
                                bji = b[ji + j];
/*<       u1 = bjf + bji >*/
                                u1 = bjf + bji;
/*<       bjg = b(jg+j) >*/
                                bjg = b[jg + j];
/*<       u2 = bjg - 0.5 * u1 >*/
                                u2 = bjg - u1 * (float).5;
/*<       u3 = c1 * ( bjf - bji ) >*/
                                u3 = c1 * (bjf - bji);
/*<       u1 = bjg + u1 >*/
                                u1 = bjg + u1;
/*<       b(jg+j) = bjc >*/
                                b[jg + j] = bjc;
/*<       a(jc+j) = t1 >*/
                                a[jc + j] = t1;
/*<       b(jc+j) = u1 >*/
                                b[jc + j] = u1;
/*<       a(jf+j) = co1*(t2-u3) - si1*(u2+t3) >*/
                                a[jf + j] = co1 * (t2 - u3) - si1 * (u2 + t3);
/*<       b(jf+j) = si1*(t2-u3) + co1*(u2+t3) >*/
                                b[jf + j] = si1 * (t2 - u3) + co1 * (u2 + t3);
/*<       a(ji+j) = co2*(t2+u3) - si2*(u2-t3) >*/
                                a[ji + j] = co2 * (t2 + u3) - si2 * (u2 - t3);
/*<       b(ji+j) = si2*(t2+u3) + co2*(u2-t3) >*/
                                b[ji + j] = si2 * (t2 + u3) + co2 * (u2 - t3);
/*<       j = j + jump >*/
                                j += *jump;
/*<   440 continue >*/
/* L440: */
                            }
/* -----(end of loop across transforms) */
/*<       ja = ja + jstepx >*/
                            ja += jstepx;
/*<       if (ja.lt.istart) ja = ja + ninc >*/
                            if (ja < istart) {
                                ja += ninc;
                            }
/*<   445 continue >*/
/* L445: */
                        }
/*<   450 continue >*/
/* L450: */
                    }
/*<   460 continue >*/
/* L460: */
                }
/* -----( end of double loop for this k ) */
/*<       kk = kk + 2*la >*/
                kk += la << 1;
/*<   470 continue >*/
/* L470: */
            }
/* -----( end of loop over values of k ) */
/*<       la = 3*la >*/
            la *= 3;
/*<   480 continue >*/
/* L480: */
        }
/* -----( end of loop on type II radix-3 passes ) */
/* -----( nvex transforms completed) */
/*<   490 continue >*/
L490:
/*<       istart = istart + nvex * jump >*/
        istart += nvex * *jump;
/*<   500 continue >*/
/* L500: */
    }
/* -----( end of loop on blocks of transforms ) */

/*<       return >*/
    return 0;
/*<       end >*/
} /* gpfa3f_ */

#ifdef __cplusplus
        }
#endif