~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/covira/ims.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /*
  2  * ims.c
  3  *
  4  * imstack infrastructure functions
  5  *
  6  */
  7 
  8 #include <tina/all_tina.h>
  9 #include <tina/brain.h>
 10 #include <tina/brainfuncs.h>
 11 
 12 static Imstack *ims = NULL; /* working image stack */
 13 
 14 /*
 15 Make an image stack with slices lz <= z < uz,
 16 Images must be of variable type vtype (uchar_v, etc).
 17 */
 18 void ims_make(Vartype vtype, int lz, int uz, double zscale)
 19 {
 20     imstack_free(ims);
 21     ims = imstack_make(vtype, lz, uz, zscale);
 22 }
 23 
 24 /*
 25 Reverse order of image stack.
 26 */
 27 void ims_reverse_image_order()
 28 {
 29     int z, lz, uz;
 30     Imrect **temp;
 31     lz = ims->lz;
 32     uz = ims->uz;
 33     temp = tvector_alloc(lz, uz, Imrect *);
 34     for (z = lz; z < uz; z++)
 35         temp[z] = ims->slice[z]->im;
 36     for (z = lz; z < uz; z++)
 37         ims->slice[z]->im = temp[uz - 1 + lz - z];
 38     tvector_free(temp, lz, Imrect *);
 39 }
 40 
 41 /*
 42 Enhance image intensities in image stack. Intensities in range low to
 43 high are gamma factored and scaled to 0-255. Outside this range they
 44 are cut-off.
 45 */
 46 void ims_rescale_images(double gamma, double low, double high)
 47 {
 48     int z;
 49     for (z = ims->lz; z < ims->uz; z++)
 50     {
 51         Imrect *im = ims->slice[z]->im;
 52         im_gamma_scale_range_inplace(im, gamma, low, high, 0.0, 255.0, 0.0, 255.0);
 53     }
 54 }
 55 
 56 /*
 57 Get global image stack.
 58 */
 59 Imstack *ims_imstack_get(void)
 60 {
 61     return ims;
 62 }
 63 
 64 /*
 65 Set global image stack.
 66 */
 67 void ims_imstack_set(Imstack * new_imstack)
 68 {
 69     imstack_free(ims);
 70     ims = new_imstack;
 71 }
 72 
 73 /*
 74 Get active VOI number.
 75 */
 76 int ims_nvoi_get()
 77 {
 78     if (ims == NULL)
 79         return(-1);
 80     return(ims->nvoi);
 81 }
 82 
 83 /*
 84 Set active VOI number (< NVOI in brain.h).
 85 */
 86 void ims_nvoi_set(int nvoi)
 87 {
 88     if (ims == NULL)
 89         return;
 90     ims->nvoi = nvoi;
 91 }
 92 
 93 /*
 94 Get active VOI in working slice.
 95 */
 96 Voi *ims_voi_get(void)
 97 {
 98     return (VOI(ims));
 99 }
100 
101 /*
102 Remove contour defined for active voi in current slice.
103 */
104 void ims_voi_empty(void)
105 {
106     if (ims != NULL)
107         voi_empty(VOI(ims));
108 }
109 
110 /*
111 Translate active VOI contour  in working slice by dp.
112 */
113 void ims_voi_shift(Vec2 dp)
114 {
115     voi_shift(VOI(ims), dp);
116 }
117 
118 /*
119 Get lower and (exclusive) upper slice numbers for
120 which current VOI has contours defined.
121 A side effect is to realise all the spline represenations
122 in the VOI.
123 */
124 Bool ims_voi_bounds(int *lz, int *uz)
125 {
126     Bool flag = false;
127     int oldz, z;
128     if (ims == NULL)
129         return(false);
130     oldz = ims_z_get();
131     /* realise all splines */
132     for (z = ims->lz; z < ims->uz; z++)
133     {
134         ims_z_set(z);
135         if(ims_spline_get() != NULL)
136             flag = true;
137     }
138     if(!flag)
139     {
140         *lz = *uz = 0;
141         return(false);
142     }
143 
144     /* find first and last defined splines */
145     for (z = ims->lz; z < ims->uz; z++)
146         if (ims->slice[z]->voi[ims->nvoi]->spline != NULL)
147         {
148             *lz = z;
149             break;
150         }
151     for (z = ims->uz-1; z >= ims->lz; z--)
152         if (ims->slice[z]->voi[ims->nvoi]->spline != NULL)
153         {
154             *uz = z+1;
155             break;
156         }
157     ims_z_set(oldz);
158     return(true);
159 }
160 
161 void ims_voi_print(FILE *fp)
162 {
163     int oldz, z;
164     Spline2 *spline;
165     if (ims == NULL)
166        return;
167     oldz = ims_z_get();
168     for (z = ims->lz; z < ims->uz; z++)
169     {
170         ims_z_set(z);
171         if((spline = ims_spline_get()) != NULL)
172         {
173            fprintf(fp,"slicez %d \n", z);
174            spline2_print(fp, spline);
175         } 
176     }
177     ims_z_set(oldz);
178 }
179 
180 void ims_voi_read(FILE *fp)
181 {
182     int oldz, z, num;
183     char temp[24];
184     Spline2 *spline;
185     if (ims == NULL)
186        return;
187     oldz = ims_z_get();
188     num = fscanf(fp,"%*s %d",&z);
189     do
190     {
191          if (num != 0)
192          {
193             ims_z_set(z);
194             spline = spline2_read(fp);
195             ims_spline_set(spline);
196             ims_spline_changed();
197             num = fscanf(fp,"%*s %d",&z);
198          }
199     } while (num>0);
200     ims_z_set(oldz);
201 }
202 
203 /*
204 Set up cursor in central (x, y, z) position.
205 */
206 void ims_init_cursor(void)
207 {
208     if (ims == NULL)
209         return;
210     ims->x = (ims->lx + ims->ux) / 2;
211     ims->y = (ims->ly + ims->uy) / 2;
212     ims->z = (ims->lz + ims->uz) / 2;
213 }
214 
215 /*
216 Get x-coordinate of 3D-cursor.
217 */
218 int ims_x_get(void)
219 {
220     if (ims == NULL)
221         return -1;
222     else
223         return ims->x;
224 }
225 
226 /*
227 Set x-coordinate of 3D-cursor.
228 */
229 void ims_x_set(int newx)
230 {
231     if (ims != NULL)
232     {
233         newx = MAX(newx, ims->lx);
234         newx = MIN(newx, ims->ux - 1);
235         ims->x = newx;
236     }
237 }
238 
239 /*
240 Get y-coordinate of 3D-cursor.
241 */
242 int ims_y_get(void)
243 {
244     if (ims == NULL)
245         return -1;
246     else
247         return ims->y;
248 }
249 
250 /*
251 Set y-coordinate of 3D-cursor.
252 */
253 void ims_y_set(int newy)
254 {
255     if (ims != NULL)
256     {
257         newy = MAX(newy, ims->ly);
258         newy = MIN(newy, ims->uy - 1);
259         ims->y = newy;
260     }
261 }
262 
263 /*
264 Get z-coordinate of 3D-cursor. This is the number of the current
265 working slice.
266 */
267 int ims_z_get(void)
268 {
269     if (ims == NULL)
270         return -1;
271     else
272         return ims->z;
273 }
274 
275 /*
276 Set z-coordinate of 3D-cursor, see ims_z_get.
277 working slice.
278 */
279 void ims_z_set(int newz)
280 {
281     if (ims != NULL)
282     {
283         newz = MAX(newz, ims->lz);
284         newz = MIN(newz, ims->uz - 1);
285         ims->z = newz;
286     }
287 }
288 
289 /*
290 Get lowest slice number in image stack (inclusive bound).
291 */
292 int ims_lz_get(void)
293 {
294     if (ims == NULL)
295         return -1;
296     else
297         return ims->lz;
298 }
299 
300 /*
301 Get highest slice number+1 in image stack (exclusive bound).
302 */
303 int ims_uz_get(void)
304 {
305     if (ims == NULL)
306         return -1;
307     else
308         return ims->uz;
309 }
310 
311 /*
312 Get inter-slice distance.
313 Scaled z-distances will only be used in explicitly 3D represenations,
314 such as graphics, spline surfaces etc., see ims_zscaled_get.
315 */
316 double ims_zscale_get(void)
317 {
318     if (ims == NULL)
319         return 1.0;
320     else
321         return ims->zscale;
322 }
323 
324 /*
325 Set inter-slice distance, see ims_zscale_get.
326 */
327 void ims_zscale_set(double zscale)
328 {
329     if (ims == NULL)
330         return;
331     ims->zscale = zscale;
332 }
333 
334 /*
335 Get slice spatial z-position (accounts for inter-slice distance factor),
336 see ims_zscale_get.
337 */
338 double ims_zscaled_get(void)
339 {
340     int z = ims_z_get();
341     return (imstack_zscaled(ims, z));
342 }
343 
344 /*
345 Set working slice one slice lower.
346 */
347 void ims_down(void)
348 {
349     if (ims != NULL)
350         ims->z = MIN(ims->z + 1, ims->uz - 1);
351 }
352 
353 /*
354 Set working slice one slice higher.
355 */
356 void ims_up(void)
357 {
358     if (ims != NULL)
359         ims->z = MAX(ims->z - 1, ims->lz);
360 }
361 
362 /*
363 Set working slice to bottom of stack.
364 */
365 void ims_bottom(void)
366 {
367     if (ims != NULL)
368         ims->z = ims->uz - 1;
369 }
370 
371 /*
372 Set working slice to top of stack.
373 */
374 void ims_top(void)
375 {
376     if (ims != NULL)
377         ims->z = ims->lz;
378 }
379 
380 /*
381 Mark current working slice for later reference.
382 */
383 void ims_voi_mark(void)
384 {
385     ims->zmark = ims->z;
386 }
387 
388 /*
389 Copy current VOI from marked slice (see ims_voi_mark) into current slice.
390 */
391 void ims_voi_import(void)
392 {
393     if (ims == NULL || ims->zmark < ims->lz || ims->zmark >= ims->uz || ims->zmark == ims->z)
394         return;
395     slice_voi_import(SLICE(ims), ims->slice[ims->zmark], ims->nvoi);
396 }
397 
398 /*
399 Get image data from current slice.
400 */
401 Imrect *ims_image_get(void)
402 {
403     if (ims == NULL)
404         return NULL;
405     else
406         return slice_image_get(SLICE(ims));
407 }
408 
409 /*
410 Set image data in current slice.
411 */
412 void ims_image_set(Imrect * im)
413 {
414     Slice *slice;
415     Imregion *roi;
416 
417     if (ims != NULL)
418     {
419         slice = SLICE(ims);
420         slice_image_set(slice, im);
421 
422         if (im != NULL)
423         {
424             roi = im->region;
425             ims->lx = MIN(ims->lx, roi->lx);
426             ims->ux = MAX(ims->ux, roi->ux);
427             ims->ly = MIN(ims->ly, roi->ly);
428             ims->uy = MAX(ims->uy, roi->uy);
429         }
430     }
431 }
432 
433 /*
434 Get snake potential image from current slice.
435 */
436 Imrect *ims_pot_get(void)
437 {
438     Imrect *pot;
439     Slice *slice;
440     double s, l, w;
441     if (ims == NULL)
442         return NULL;
443     slice = SLICE(ims);
444     pot = slice_pot_get(slice);
445     if(pot != NULL || slice->im == NULL)
446         return(pot);
447     s = ims_scale_get();
448     l = ims_level_get();
449     w = ims_window_get();
450     pot = canny_window_pot(slice->im, s, l, w);
451     slice_pot_set(slice, pot);
452     return(pot);
453 }
454 
455 /*
456 Get snake potential image from current slice.
457 */
458 void ims_pot_set(Imrect * pot)
459 {
460     if (ims == NULL)
461         return;
462     slice_pot_set(SLICE(ims), pot);
463 }
464 
465 /*
466 Free snake potentials in all files.
467 */
468 void ims_pots_free()
469 {
470     int z, lz, uz, oldz;
471     if(ims == NULL)
472         return;
473     oldz = ims_z_get();
474     lz = ims_lz_get();
475     uz = ims_uz_get();
476     for(z = lz; z < uz; z++)
477     {
478         ims_z_set(z);
479         ims_pot_set(NULL);
480     }
481     ims_z_set(oldz);
482 }
483         
484 /*
485 Get 2D point string representing active VOI in current slice.
486 */
487 Tstring *ims_string_get(void)
488 {
489     if (ims == NULL)
490         return NULL;
491     else
492         return voi_string_get(VOI(ims));
493 }
494 
495 /*
496 Set 2D point string representing active VOI in current slice.
497 */
498 void ims_string_set(Tstring * string)
499 {
500     if (ims != NULL)
501         voi_string_set(VOI(ims), string);
502 }
503 
504 /*
505 Notify change in string for active VOI in current slice.
506 Other representations will then be recomputed on ims_get_... request.
507 */
508 void ims_string_changed(void)
509 {
510     if (ims != NULL)
511         voi_string_changed(VOI(ims));
512 }
513 
514 /*
515 Get B-spline representing active VOI in current slice.
516 If not available spline will be computed by recursive fit
517 from string or snake.
518 */
519 Spline2 *ims_spline_get(void)
520 {
521     if (ims == NULL)
522         return NULL;
523     else
524         return voi_spline_get(VOI(ims));
525 }
526 
527 /*
528 Get B-spline representing active VOI in current slice. Forces
529 fitted spline to have n control points.
530 */
531 Spline2 *ims_fixed_spline_get(int n)
532 {
533     if (ims == NULL)
534         return NULL;
535     else
536         return voi_fixed_spline_get(VOI(ims), n);
537 }
538 
539 /*
540 Set B-spline representing active VOI in current slice.
541 */
542 void ims_spline_set(Spline2 * spline)
543 {
544     if (ims != NULL)
545         voi_spline_set(VOI(ims), spline);
546 }
547 
548 /*
549 Notify change in spline for active VOI in current slice.
550 Other representations will then be recomputed on ims_get_... request.
551 */
552 void ims_spline_changed(void)
553 {
554     if (ims != NULL)
555         voi_spline_changed(VOI(ims));
556 }
557 
558 /*
559 Get snake representing active VOI in current slice.
560 If not available snake will be computed from string or snake.
561 Knots will be chosen to ensure average angle at vertices is less than
562 asnake.
563 */
564 Snake *ims_snake_get(void)
565 {
566     if (ims == NULL)
567         return NULL;
568     else
569         return voi_snake_get(VOI(ims));
570 }
571 
572 /*
573 Set snake representing active VOI in current slice.
574 */
575 void ims_snake_set(Snake * snake)
576 {
577     if (ims != NULL)
578         voi_snake_set(VOI(ims), snake);
579 }
580 
581 /*
582 Notify change in snake for active VOI in current slice.
583 Other representations will then be recomputed on ims_get_... request.
584 */
585 void ims_snake_changed(void)
586 {
587     if (ims != NULL)
588         voi_snake_changed(VOI(ims));
589 }
590 
591 /*
592 Get region of interest of image in current working slice.
593 */
594 Imregion *ims_roi_get(void)
595 {
596     Imregion *roi;
597 
598     if (ims == NULL)
599         return NULL;
600     else
601     {
602         roi = roi_alloc(ims->lx, ims->ly, ims->ux, ims->uy);
603         return roi;
604     }
605 }
606 
607 /*
608 Return transverse image in (y, z) plane at current x-cursor position.
609 */
610 Imrect *ims_xslice(void)
611 {
612     if (ims == NULL)
613         return NULL;
614     else
615         return imstack_xslice(ims, ims->x);
616 }
617 
618 /*
619 Return transverse image in (x, z) plane at current y-cursor position.
620 */
621 Imrect *ims_yslice(void)
622 {
623     if (ims == NULL)
624         return NULL;
625     else
626         return imstack_yslice(ims, ims->y);
627 }
628 
629 /*
630 Return current slice. Note difference from ims_xslice and ims_yslice.
631 */
632 Imrect *ims_zslice(void)
633 {
634     if (ims == NULL)
635         return NULL;
636     else
637         return SLICE(ims)->im;
638 }
639 
640 /*
641 Set force law to be used to calculate forces at snake knots.
642 Defaults is snake_gradient_forces (see snake_run.c).
643 */
644 void ims_snake_forcelaw_set(void (*forcelaw)())
645 {
646     int i, z;
647     if (ims == NULL)
648         return;
649     for (z = ims->lz; z < ims->uz; z++)
650     {
651         Slice *slice = ims->slice[z];
652         for (i = 0; i < NVOI; i++)
653             if (slice->voi[i]->snake != NULL)
654                 slice->voi[i]->snake->forcelaw = forcelaw;
655     }
656 }
657 
658 /*
659 Compute potential image in each slice.
660 */
661 void ims_pots_get(void)
662 {
663     int z, zold;
664     if(ims == NULL)
665         return;
666     zold = ims_z_get();
667     for (z = ims->lz; z < ims->uz; z++)
668     {
669         ims_z_set(z);
670         ims_pot_get();
671     }
672     ims_z_set(zold);
673 }
674 
675 /*
676 Run current snake in translation only mode in potential image
677 for given number of steps.
678 */
679 double ims_snake_run_trans(int steps)
680 {
681     double d;
682     Snake *snake = ims_snake_get();
683     Imrect *im = ims_pot_get();
684 
685     if (ims == NULL || snake==NULL || im==NULL) /* NAT 27/2/2000 */
686         return 0.0;
687     else
688     {
689         d = snake_run_trans(snake, im, steps);
690         ims_snake_changed();
691         return d;
692     }
693 }
694 
695 /*
696 Run current snake in rotation only mode in potential image
697 for given number of steps.
698 */
699 double ims_snake_run_rot(int steps)
700 {
701     double d;
702     Snake *snake = ims_snake_get();
703     Imrect *im = ims_pot_get();
704 
705     if (ims == NULL || snake==NULL || im==NULL) /* NAT 27/2/2000 */
706         return 0.0;
707     else
708     {
709         d = snake_run_rot(snake, im, steps);
710         ims_snake_changed();
711         return d;
712     }
713 }
714 
715 /*
716 Run current snake in scale only mode in potential image
717 for given number of steps.
718 */
719 double ims_snake_run_scale(int steps)
720 {
721     double d;
722     Snake *snake = ims_snake_get();
723     Imrect *im = ims_pot_get();
724 
725     if (ims == NULL || snake==NULL || im==NULL) /* NAT 27/2/2000 */
726         return 0.0;
727     else
728     {
729         d = snake_run_scale(snake, im, steps);
730         ims_snake_changed();
731         return d;
732     }
733 }
734 
735 /*
736 Run current snake in rotation, translation and scale only mode in
737 potential image for given number of steps.
738 */
739 double ims_snake_run_rts(int steps)
740 {
741     double d;
742     Snake *snake = ims_snake_get();
743     Imrect *im = ims_pot_get();
744 
745     if (ims == NULL || snake==NULL || im==NULL) /* NAT 27/2/2000 */
746         return 0.0;
747     else
748     {
749         d = snake_run_rts(snake, im, steps);
750         ims_snake_changed();
751         return d;
752     }
753 }
754 
755 /*
756 Run current snake in all-deformations-allowed mode in
757 potential image for given number of steps.
758 */
759 double ims_snake_run_all(double alpha, double beta, int steps)
760 {
761     double d;
762     Snake *snake = ims_snake_get();
763     Imrect *im = ims_pot_get();
764 
765     if (ims == NULL || snake==NULL || im==NULL) /* NAT 27/2/2000 */
766         return 0.0;
767     else
768     {
769         d = snake_run_all(snake, im, alpha, beta, steps);
770         ims_snake_changed();
771         return d;
772     }
773 }
774 
775 /*
776 Get intensity data from image and put it in stats.
777 data set up by calling function, used by apply_...
778 functions.
779 */
780 static void stats_scan(int y, int lx, int ux, void **data)
781 {
782     Stats *stats = data[0];
783     Imrect *im = data[1];
784     int x, pix;
785     for (x = lx; x <= ux; x++)
786     {
787         stats->n++;
788         pix = im_get_pix(im, y, x);
789         stats->tot1 += pix;
790         stats->tot2 += pix * pix;
791     }
792 }
793 
794 /*
795 Get statistics (mean and standard deviation) of intensity image
796 in interior of current snake, and store them for
797 use by statistical snake.
798 */
799 void ims_stat_snake_init()
800 {
801     int n;
802     Snake *snake = ims_snake_get();
803     Imrect *im = ims_image_get();
804     Stats *stats;
805     void *data[2];
806 
807     if (!snake || !im)
808         return;
809     stats = stats_make();
810     data[0] = stats;
811     data[1] = im;
812     apply_inside_poly(snake->n, snake->p, im->width, im->height, stats_scan, data);
813     if ((n = stats->n) > 1)
814     {
815         stats->mean = stats->tot1 / n;
816         stats->var = (stats->tot2 - sqr(stats->tot1) / n) / (n - 1);
817         stats->sd = sqrt(stats->var);
818     }
819     ims_stats_set(stats);
820     printf("N: %d\tMean: %.3f\tSD: %.3f\n", stats->n, stats->mean, stats->sd);
821 }
822 
823 /*
824 Run current snake as statistical region grower,
825 in all-deformations-allowed mode,
826 for given number of steps.
827 Intensity is accepted if it is within k_s_d standard deviations of mean.
828 If flag is true, snake will dynamically insert and delete knots, keeping
829 inter-knot gap between min_gap and max_gap.
830 The knot which moves furthest will move a distance dt.
831 */
832 void ims_stat_snake_run(double alpha, double beta, double pressure,
833           double k_s_d, int steps, int flag,
834           int min_gap, int max_gap, double dt)
835 {
836     Stats *stats = ims_stats_get();
837     Imrect *im = ims_image_get(), *accum;
838     Tstring *ss = ss_from_snake(ims_snake_get());
839     int i;
840 
841     if (ss == NULL || im == NULL)
842         return;
843 
844     accum = accum_make(im, ss);
845     if(flag)
846     {
847         stat_snake_mingap(ss, accum, min_gap);
848         stat_snake_maxgap(ss, accum, min_gap, max_gap);
849     }
850     else
851         stat_snake_mingap(ss, accum, 1);
852     for (i = 0; i < steps; i++)
853     {
854         stat_snake_step(ss, accum, im, stats, alpha, beta, pressure, k_s_d, dt);
855         if(flag)
856         {
857             stat_snake_mingap(ss, accum, min_gap);
858             stat_snake_maxgap(ss, accum, min_gap, max_gap);
859         }
860     }
861     ims_snake_set(ss_to_snake(ss)); ss_free(ss);
862     ims_snake_changed();
863     im_free(accum);
864 }
865 
866 /*
867 Run current snake as statistical region grower,
868 in tanslation only mode, for given number of steps.
869 Intensity is accepted if it is within k_s_d standard deviations of mean.
870 If flag is true, snake will dynamically insert and delete knots, keeping
871 inter-knot gap between min_gap and max_gap.
872 The knot which moves furthest will move a distance dt.
873 */
874 void ims_stat_snake_trans(double alpha, double beta, double pressure,
875           double k_s_d, int steps, double dt)
876 {
877     Stats *stats = ims_stats_get();
878     Imrect *im = ims_image_get();
879     Tstring *ss = ss_from_snake(ims_snake_get());
880     int i;
881     if (ss == NULL || im == NULL)
882         return;
883     for (i = 0; i < steps; i++)
884         stat_snake_trans(ss, im, stats, alpha, beta, pressure, k_s_d, dt);
885     ims_snake_set(ss_to_snake(ss)); ss_free(ss);
886     ims_snake_changed();
887 }
888 
889 /*
890 Use data in current VOI to build a 3D B-spline surface
891 model. This model will be a tube (closed in (x,y) planes, open in z).
892 */
893 Splinesurf3 *ims_splinesurf3_make()
894 {
895     int z, zold, lz, uz, u, v, nu, nv;
896     void *surf;
897     Spline2 *spline;
898     Vec3 **data;
899 
900     if (ims == NULL)
901         return (NULL);
902     imstack_uniform_splines(ims);
903     zold = ims_z_get();
904 
905     spline = NULL;
906     for (z = ims->lz; z < ims->uz && spline == NULL; z++)
907     {
908         ims_z_set(z);
909         spline = ims_spline_get();
910     }
911     lz = z - 1;
912     spline = NULL;
913     for (z = ims->uz - 1; z >= ims->lz && spline == NULL; z--)
914     {
915         ims_z_set(z);
916         spline = ims_spline_get();
917     }
918     uz = z + 2;
919 
920     if (lz >= uz)
921         return (NULL);
922     ims_z_set(lz);
923     spline = ims_spline_get();
924     nu = spline->n;
925     nv = uz - lz;
926     surf = splinesurf3_make(SPLINE_PERIODIC, SPLINE_NATURAL, nu, nv);
927     data = tarray_alloc(0, 0, nu, nv, Vec3);
928 
929     for (z = lz; z < uz; z++)
930     {
931         v = z-lz;
932         ims_z_set(z);
933         spline = ims_spline_get();
934         for (u = 0; u < nu; u++)
935         {
936             Vec2 p = spline2_eval(spline, u);
937             data[u][v] = vec3(vec3_x(p), vec3_y(p), ims_zscaled_get());
938         }
939     }
940     splinesurf3_interpolate(surf, data);
941     tarray_free(data, 0, 0, nu, nv, Vec3);
942 
943     ims_z_set(zold);
944     return (surf);
945 }
946 
947 /*
948 Use all data in current VOI to interpolate a contour
949 into current slice.
950 */
951 void ims_interpolate_spline(void)
952 {
953     double z;
954  
955     if (ims != NULL)
956     {
957         Spline2 *spline;
958         z = ims_z_get();
959         imstack_uniform_splines(ims);
960         spline = imstack_interpolate_spline(ims, z);
961         ims_spline_set(spline);
962     }
963 }
964  
965 /*
966 Use all data in current VOI between lz and (exclusive) uz
967 to interpolate contours into allslice between these values.
968 */
969 Bool ims_interpolate_all_splines_between(int lz, int uz)
970 {
971     Spline2 **spline;
972     int oldz, z;
973 
974     if (ims == NULL)
975         return(false);
976     imstack_uniform_splines(ims);
977     if(lz >= uz)
978         return(false);
979 
980     oldz = ims_z_get();
981 
982     /* interpolate all null internal splines */
983     spline = tvector_alloc(lz, uz, Spline2 *);
984     for (z = lz; z < uz; z++)
985     {
986         ims_z_set(z);
987         if (ims_spline_get() == NULL)
988             spline[z] = imstack_interpolate_spline(ims, (double) z);
989     }
990 
991     /** after interpolation save splines **/
992     for (z = lz; z < uz; z++)
993     {
994         ims_z_set(z);
995         if (ims_spline_get() == NULL)
996         {
997             ims_spline_set(spline[z]);
998             ims_spline_changed();
999         }
1000     }
1001 
1002     tvector_free(spline, lz, Spline2 *);
1003 
1004     ims_z_set(oldz);
1005     return(true);
1006 }
1007 
1008 /*
1009 Free image stack and all associated memory.
1010 */
1011 void ims_free()
1012 {
1013     imstack_free(ims);
1014 }
1015 
1016 /*
1017 Gets the tube (transverse boundary) structure to be used
1018 for tube extrapolation.
1019 */
1020 Tube *ims_tube_get()
1021 {
1022     if(ims == NULL)
1023         return(NULL);
1024     return(ims->tube);
1025 }
1026 
1027 /*
1028 Sets the tube (transverse boundary) structure to be used
1029 for tube extrapolation.
1030 */
1031 void ims_tube_set(Tube *newtube)
1032 {
1033     if(ims == NULL)
1034         return;
1035     tube_free(ims->tube);
1036     ims->tube = newtube;
1037 }
1038 
1039 /*
1040 Sets the tube (transverse boundary) structure to be used
1041 for tube extrapolation, using data from a 3D point string.
1042 This string must be closed (LOOP type) and contain points
1043 in all intermediate slices. If the top and bottom slice
1044 contain only one point they are ignored.
1045 */
1046 void ims_string_tube_set(Tstring *str)
1047 {
1048     List *a1, *a2, *b1, *b2, *c1, *c2, *ptr;
1049     int lz, uz, z1, z2;
1050     Tube *tube;
1051 
1052     if (str == NULL)
1053         return;
1054     lz = ims_lz_get();
1055     uz = ims_uz_get();
1056     z1 = MAXINT;
1057     z2 = -MAXINT;
1058     for(ptr = str->start;; ptr = ptr->next)
1059     {
1060         int z = vec3_z(*(Vec3 *)ptr->to);
1061         if(z1 > z)
1062         {
1063             a1 = ptr;
1064             z1 = z;
1065         }
1066         if(z2 < z)
1067         {
1068             a2 = ptr;
1069             z2 = z;
1070         }
1071         if(ptr == str->end)
1072             break;
1073     }
1074     z1 = MAX(z1, lz);
1075     z2 = MIN(z2, uz-1);
1076     if(z1 == z2)
1077     {
1078         ims_tube_set(NULL);
1079         return;
1080     }
1081 
1082     for(ptr = a1; ;ptr = ptr->next)
1083     {
1084         int z = vec3_z(*(Vec3 *)ptr->to);
1085         if(z > z1)
1086         {
1087             b1 = ptr->last;
1088             break;
1089         }
1090     }
1091     for(ptr = a1; ;ptr = ptr->last)
1092     {
1093         int z = vec3_z(*(Vec3 *)ptr->to);
1094         if(z > z1)
1095         {
1096             c1 = ptr->next;
1097             break;
1098         }
1099     }
1100     for(ptr = a2; ;ptr = ptr->last)
1101     {
1102         int z = vec3_z(*(Vec3 *)ptr->to);
1103         if(z < z2)
1104         {
1105             b2 = ptr->next;
1106             break;
1107         }
1108     }
1109     for(ptr = a2; ;ptr = ptr->next)
1110     {
1111         int z = vec3_z(*(Vec3 *)ptr->to);
1112         if(z < z2)
1113         {
1114             c2 = ptr->last;
1115             break;
1116         }
1117     }
1118     if(b1 == c1)
1119     {
1120         b1 = b1->next;
1121         c1 = c1->last;
1122         z1 = vec3_z(*(Vec3 *)b1->to);
1123     }
1124     if(b2 == c2)
1125     {
1126         b2 = b2->last;
1127         c2 = c2->next;
1128         z2 = vec3_z(*(Vec3 *)b2->to);
1129     }
1130     tube = tube_make(z1, z2+1);
1131     for(ptr = b1; ; ptr = ptr->next)
1132     {
1133         Vec3 p = *(Vec3 *)ptr->to;
1134         int z = vec3_z(p);
1135         tube->p1[z] = vec2(vec3_x(p), vec3_y(p));
1136         if(ptr == b2)
1137             break;
1138     }
1139     for(ptr = c1; ; ptr = ptr->last)
1140     {
1141         Vec3 p = *(Vec3 *)ptr->to;
1142         int z = vec3_z(p);
1143         tube->p2[z] = vec2(vec3_x(p), vec3_y(p));
1144         if(ptr == c2)
1145             break;
1146     }
1147     ims_tube_set(tube);
1148 }
1149 
1150 /*
1151 Define a tube transverse boundary using 2 B-splines
1152 in a given sagittal plane.
1153 */
1154 void ims_sagit_tube_set(Spline2 * spline1, Spline2 *spline2, double x)
1155 {
1156     Vec2 p1, p2;
1157     int z, lz, uz, z1, z2;
1158     Spline *zspline;
1159     Tube *tube;
1160 
1161     if (spline1 == NULL || spline2 == NULL)
1162         return;
1163     lz = ims_lz_get();
1164     uz = ims_uz_get();
1165 
1166     p1 = spline2_eval(spline1, 0);
1167     p2 = spline2_eval(spline1, spline1->n - 1.0);
1168     z1 = ROUND(vec2_y(p1));
1169     z2 = ROUND(vec2_y(p2)) + 1;
1170     if (z1 > z2)
1171         SWAP(int, z1, z2);
1172     if (z1 < lz)
1173         z1 = lz;
1174     if (z2 > uz)
1175         z2 = uz;
1176 
1177     lz = z1; uz = z2;
1178     p1 = spline2_eval(spline2, 0);
1179     p2 = spline2_eval(spline2, spline2->n - 1.0);
1180     z1 = ROUND(vec2_y(p1));
1181     z2 = ROUND(vec2_y(p2)) + 1;
1182     if (z1 > z2)
1183         SWAP(int, z1, z2);
1184     if (z1 < lz)
1185         z1 = lz;
1186     if (z2 > uz)
1187         z2 = uz;
1188 
1189     tube = tube_make(z1, z2);
1190 
1191     zspline = spline1->y;
1192     for (z = z1; z < z2; z++)
1193     {
1194         double t;
1195         Vec2 p, q;
1196         t = spline_param(zspline, z);
1197         p = spline2_eval(spline1, t);
1198         q = vec2(x, vec2_x(p));
1199         tube->p1[z] = q;
1200     }
1201     zspline = spline2->y;
1202     for (z = z1; z < z2; z++)
1203     {
1204         double t;
1205         Vec2 p, q;
1206         t = spline_param(zspline, z);
1207         p = spline2_eval(spline2, t);
1208         q = vec2(x, vec2_x(p));
1209         tube->p2[z] = q;
1210     }
1211     ims_tube_set(tube);
1212 }
1213 
1214 /*
1215 Find rotation matrix, translation vector, and scaling,
1216 that takes line (p1, p2) to line (q1, q2)
1217 q = srp+t
1218 */
1219 void rts_get(Vec2 p1, Vec2 p2, Vec2 q1, Vec2 q2,
1220     Mat2 *r, Vec2 *t, double *s)
1221 {
1222     double theta = vec2_angle(vec2_diff(q2, q1), vec2_diff(p2, p1));
1223     *r = rot2(theta);
1224     *s = vec2_dist(q1, q2)/vec2_dist(p1, p2);
1225     *t = vec2_diff(q1, vec2_times(*s, mat2_vprod(*r, p1)));
1226 }
1227  
1228 /*
1229 Use current tube transverse boundary representation to
1230 extrapolate all current contours in current VOI
1231 into all slices for which tube data exists.
1232 The original contour is translated, rotated, and scaled to match
1233 the boundary data.
1234 */
1235 
1236 void ims_tube_interp()
1237 {
1238     Tube *tube;
1239     Spline2 *spline;
1240     Mat2 r;
1241     Vec2 p1, p2, q1, q2, t;
1242     double s;
1243     int z, lz, uz, z1, z2, zold = ims_z_get();
1244  
1245     ims_voi_bounds(&lz, &uz);
1246     ims_interpolate_all_splines_between(lz, uz);
1247 
1248     tube = ims_tube_get();
1249     if(tube == NULL)
1250         return;
1251     z1 = tube->z1;
1252     z2 = tube->z2;
1253  
1254     for(z = lz; z < uz; z++)
1255     {
1256         ims_z_set(z);
1257         q1 = tube->p1[z];
1258         q2 = tube->p2[z];
1259         spline = ims_spline_get();
1260         p1 = spline2_closest(spline, q1);
1261         p2 = spline2_closest(spline, q2);
1262         rts_get(p1, p2, q1, q2, &r, &t, &s);
1263         spline2_rts(spline, r, t, s);
1264         ims_spline_changed();
1265     }
1266 
1267     ims_z_set(lz); spline = ims_spline_get();
1268     for(z = lz-1; z >= z1; z--)
1269     {
1270         ims_z_set(z+1); spline = spline2_copy(ims_spline_get());
1271         ims_z_set(z);
1272         ims_z_set(z);
1273         q1 = tube->p1[z];
1274         q2 = tube->p2[z];
1275         p1 = spline2_closest(spline, q1);
1276         p2 = spline2_closest(spline, q2);
1277         rts_get(p1, p2, q1, q2, &r, &t, &s);
1278         spline2_rts(spline, r, t, s);
1279         ims_spline_set(spline2_copy(spline));
1280         ims_spline_changed();
1281     }
1282 
1283     for(z = uz-1; z < z2; z++)
1284     {
1285         ims_z_set(z-1); spline = spline2_copy(ims_spline_get());
1286         ims_z_set(z);
1287         q1 = tube->p1[z];
1288         q2 = tube->p2[z];
1289         p1 = spline2_closest(spline, q1);
1290         p2 = spline2_closest(spline, q2);
1291         rts_get(p1, p2, q1, q2, &r, &t, &s);
1292         spline2_rts(spline, r, t, s);
1293         ims_spline_set(spline);
1294         ims_spline_changed();
1295     }
1296 
1297     ims_z_set(zold);
1298 }
1299 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.